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ABSTRACT 

Two  approximate  parametric  interval  estimation  methods  for  system  reliability  using 
component  test  data  are  developed  and  evaluated.  One  method  can  be  applied  to  any 
coherent  system  with  components  which  have  exponential  failure  times  with  possibly 
different  failure  rates  and  different  mission  operating  times.  This  method  estimates  the 
ratios  of  component  failure  rates  which  are  then  used  to  develop  the  approximate  lower 
confidence  limit.  These  ratio  estimates  are  developed  with  and  without  jacknife  methods 
and  the  two  results  are  compared.  This  procedure  is  very  accurate  and  simple  to  com- 
pute, requiring  the  use  of  standard  chi-square  tables.  This  ratio  method  is  subsequently 
extended  to  coherent  systems  with  components  whose  failure  times  have  a  Weibull  dis- 
tribution. A  nearly  exact  parametric  lower  confidence  limit  for  P(X>  x)  is  developed  and 
evaluated  where  x  is  given  and  X  has  a  normal  distribution  with  unknown  mean  and 
variance.  This  procedure  is  also  simple  to  evaluate  and  requires  the  use  of  Student  t  ta- 
bles. 
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I.     BACKGROUND 

Parametric  confidence  interval  procedures  for  the  reliability  of  mechanical  systems 
are  much  less  developed  than  procedures  for  electrical  systems.  This  is  due  to  the  more 
complicated  failure  distributions  used  to  model  mechanical  hardware.  The  failure  rate 
of  a  series  system  of  independent  components,  with  exponentially  distributed  failure 
times  and  equal  mission  times,  is  the  sum  of  the  failure  component  failure  rates.  This 
property  has  permitted  development  of  numerous  methods  for  system  reliability  of  series 
and  other  coherent  systems  using  component  failure  data.  The  IVeibull  and  extreme 
value  distributions  have  been  used  in  life  testing  methods  for  both  electrical  and  me- 
chanical devices.  Several  methods  have  been  used  for  obtaining  point  estimates  of  pa- 
rameters for  these  distributions  and  for  the  reliability  function  itself.  Harter  and  Moore 
[Ref.  1:  pp.  889-901],  Mann  [Ref  2:  pp.  231-256], 

The  derivation  of  simple  confidence  limits  for  the  reliability  function  for  the  extreme 
value  distribution  with  parameters  z  and  3  has  posed  problems,  because  methods  based 
on  a  pivotal  quantity  such  as  |-  —  z)jd  are  inadequate.  Lawless  [Ref.  3:  pp.  355-364], 
Johns  and  Lieberman  [Ref.  4:  pp.  135-175]  and  Thoman,  Bain  and  Antle  [Ref.  5:  pp. 
363-372]  have  developed  nearly  exact  procedures  for  confidence  limits  for  the  reliability 
function  of  the  IVeibull  and  extreme  value  distributions.  Schneider  and  Weissfeld  [Ref. 
6:  pp.  179-186]  have  developed  interval  estimation  methods  for  percentiles  of  the  IVeibull 
and  extreme  value  distribution  based  on  censored  data.  Although  extensive  methods 
have  been  developed  for  interval  estimates  of  the  reliability  of  a  single  component  with 
the  IVeibull  or  extreme  value  failure  distribution,  very  few  parametric  interval  methods 
have  been  developed  for  system  reliability  using  component  test  data  with  IVeibull  failure 
distributions. 

Two  approximate  parametric  interval  estimation  methods  for  reliability  of  coherent 
systems  using  component  test  data  are  developed  and  evaluated  in  this  thesis.  Evalu- 
ations of  these  procedures  are  performed  using  computer  simulation  for  series  systems 
only.  One  of  these  methods  is  developed  for  the  reliability  of  a  series  system  whose 
components  have  exponential  failure  distributions  and  different  mission  times.  This 
procedure  was  found  to  be  quite  accurate  and  can  be  applied  to  coherent  systems  in 
general.  This  first  procedure  is  then  extended  to  the  case  where  components  of  the  sys- 


tern  have  Weibull  failure  distributions.  The  method  used  is  an  extension  of  a  non- 
parametric  method  developed  by  Vlyhre,  J.,  Rosenfeld,  A.  and  Saunders,  S.  [Ref.  7: 
PP.  2i: 3-227]. 

The  normal  distribution  is  used  extensively  in  some  mechanical  reliabilty  models. 
Maximum  likelihood  and  minimum  variance  unbiased  estimators  for  P(X>x0)  when 
both  (x  and  a2  are  unknown  were  developed  over  thirty  years  ago,  Lieberman  and 
Resnikoff  [Ref.  8:  pp.  457-516],  Folks  and  others  [Ref.  9:  pp.  43-50],  and  Barton  [Ref. 
10:  pp.  227-229].  Exact  interval  estimation  procedures  for  P(X>x0)  were  developed  by 
Owen  and  Hua  [Ref.  1 1:  pp. 285-31 1]  using  the  non-central  /  distribution.  Letting  X  de- 
note strength  and  Y  denote  stress,  mechanical  reliability  is  sometimes  modeled  as  the 
value  for  P(X>  Y).  Approximate  interval  estimation  procedures  for  this  probability 
when  X  and  Y  are  assumed  to  be  normal  have  been  developed  by  Church  and  Harris 
[Ref.  12:  pp.  49-54]  when  the  mean  and  variance  of  }'  are  known.  Downton  [Ref.  13:  pp. 
551-55S]  modifies  their  procedure  slightly  to  get  more  accurate  bounds  and  suggests  an 
approximate  procedure  when  the  means  and  variances  of  both  X  and  }'  are  unknown. 
Lower  confidence  intervals  for  P(X>  Y)  obtained  under  the  assumption  of  normality  for 
X  and  Y  can  lead  to  serious  error  when  either  X  or  Y  or  both  are  truncated  well  into  the 
tails.  Consequently,  P(X>  x„)  may  be  a  more  reasonable  model  of  mechanical  reliability 
where  x0  is  chosen  conservatively.  A  very  accurate  approximate  lower  confidence  limit 
procedure  for  P(X>xg)  is  developed  and  evaluated  in  this  thesis.  It  can  be  computed 
easily. 


II.     INTERVAL  ESTIMATION  PROCEDURE  -  EXPONENTIAL  CASE 

A  system  of  independent  components  is  coherent  if  an  increase  in  reliability  of  any 
one  of  its  components  does  not  cause  a  degredation  in  system  reliability.  Suppose  a 
coherent  system  has  k  components.  We  assume  that  the  failure  distribution  of  compo- 
nent /  is  exponential  with  failure  rate  >.,  .  Then  system  reliability  Rs  can  be  written  as  a 
function  of /;„ /,,     i=  1,  2,  ••  • ,  k;      i.e., 

rts(')  =  £(>-i,  >-2.  •  •  •  >'-k<  h>hi'  •  •  >h) 
where   tt  =  i,(t)  is  the  operating  time  for  component  /,  /'  =  1,  2,  •••  ,k  . 

Let  /„  =  max{i1,  /.2,  •  •  • ,  /*}  and  r,  =  >.,//m    /  =  1,  2,  •  •  • ,  k.  Then  one  can  write 

rts(')  =  gV-mr\>  •••  ,rftf  //,•••.  //J 

If  the  r,  were  known  and  /m,r(a)  were  an  upper  100(1  -  a)%  confidence  limit  for  /,„  ,  the 
corresponding  lower  confidence  limit  for  Rs(t)  would  be 

Ks(0z.(*)  =  g(>-m,U(*)>  rr  *  '  '  <''/,-  '],«••,  '/<) 

Specifically,  if  the  system  is  a  series  system  of  independent  components,  so  that 

f     V      ]  f         V      1 

K5(/)  =  exp<  -  ^  Aj/i  >  =  exp^  -  /.^  r^  ^  (2.1) 


then. 


*s(0k«)  =  exP]  -im,im£j  w 


If  h,  items  of  component  /  are  tested  until  failure,  T,  denotes  the  total  test  time  accumu- 

lated  by  all  n,  items   and  n  =  Yn,,  then    2/mX/,T,    is  y\n  .   See  Bain  and  Engelhardt  [Ref. 

1=1  i=i 

14]. 

An  upper  confidence  limit  for  /m  is 


2 
*  y-x,2n 

Am,U(x)  ~         k  ^A-^> 


i  *  i 


where  y\_n  is  the  100(1  —  a)th  percentile  of  the  xl  distribution.    Corresponding  equations 
for  truncated  testing  are  similar. 

If  the  r,  are  unknown,  the  following  methods  estimate  the  values  for  r,  from  the  data. 
One  method  uses  the  likelihood  ratio  estimate  for  r,.  The  second  method  uses  a  jacknifed 
version  of  the  first  method.  The  two  resulting  confidence  limits  RSM%)  and  Rs,u*)  >  w^^ 
and  without  jacknifing  r,  respectively,  are  compared  for  relative  accuracy.  Quenouille 
[Ref.  15:  pp.  353-360]  first  reported  a  method  for  estimating  ratios  that  reduced  the  bias 
without  increasing  the  variance.  Miller  [Ref.  16:  pp.  1-15]  gives  an  excellent  review  of  the 
jacknile  method  which  includes  a  discussion  on  the  application  of  jacknifing  to  estimat- 
ing ratios. 

A.     LOWER  CONFIDENCE  LIMIT  Rsl  WITHOUT  JACKNIFING 

In  this  case,  the  maximum  likelihood  estimate  of  the  ratio  r,  =  — -  is 


M 


(2.3) 


where  Xi  =  nijTl  and  ).m  =    max  (/.lt  •  •  • ,  Ak) .  The  resulting  approximate  upper  confi- 
dence limit  for  }.m  is 

2 
0  '?--<  In 


Z^ 


2 


k 
where      n  =  Tn,  ,  and 


n 

7=1 


Tt-yT,,      i   =    1,2,...,  A  (2.5) 


and  7  denotes  the  failure  time  of  thejth  test  for  component  i. 
The  resulting  approximate  confidence  bound  RSMa)  is  given  by 


f        <  Va    1 

B.     LOWER  CONFIDENCE  LIMIT  RSL  WITH  JACKNIFING 

The  definitions  for  /,  ,  r,  ,  n,  ,  and  7*  in  Section  A  are  also  used  in  this  section.    Let 
/,,(.)  denote  the  estimate  for  >.,  by  removing  TtJ  from  the  data  ;  i.e., 

^(.)==_^r-         J   =    1,2,  •..,«,-,      l+j  (2.7) 

and 


j    =    l,2....,nm,       l*j  (2.S) 


Then  the  jacknifed  ratio  estimate  r*  is  given  by 


(/'(•) 


A*  *A  /=! 


where        «"  =  min(«i,  n2,  ••  • ,  nk)        and 


Now  define    /.'lM  by 


2 


1^ 


/=! 


The  corresponding  confidence  bound  Rs,iM  IS  given  by 


(2.9) 


ty(.)  =  ~ J    =    1,2, ...,«  (2.10) 

Amj(>) 


2 
A*  ^i  In 

im—r —  •  <2"> 


^X(,)  =  eXP)   -/L(,)Z^'''/'(  <112) 

When  the  nit  i  =  1,2,  •»*,k  difTer  considerably,  this  jacknife  estimation  proce- 
dure can  be  unbalanced.  That  is,  the  number  of  data  points  used  to  compute  the  jacknife 
estimate  will  differ  from  one  component  to  another.  It  was  decided  to  use  this  jacknifing 
procedure  and  determine  the  effect  of  differing  sample  sizes  by  examining  the  results  of 
the  simulations. 

In  equation  (2.9),  the  jacknife  estimate  is  constructed  by  using  only  the  first  n  ob- 
servations from  each  component  to  obtain  rIJ(.}  where  n  =  mm{nlf  n2,  •  •  • ,  nk) .  Of 
course  it  is  rather  arbitrary  to  take  the  first  ri  observations  from  n,  .  since  any  n  of  the 
collection  of  n,  values  could  be  used. 

C.     SIMULATIONS  AND  RESULTS  TOR  EXPONENTIAL  CASE 

1.     Simulation 

a.  Simulation  language  and  package 

The  programming  language  used  to  simulate  this  problem  was  VS 
FORTRAN  on  an  IBM  3033.  Also  LEXPN  in  LLRAXII  was  used  to  generate  observed 
exponential  random  variatcs.    SHSORT  was  used  to  perform  the  sort  routines. 

b.  Cases  and  Input  parameters 

The  six  input  parameters  below  determine  the  conditions  for  exercising  the 
simulation  runs.  System  reliability  is  determined  by  five  parameters.  Values  of  these 
parameters  are  given  in  the  tables  that  show  the  results.  The  test  plan  simulated  was  to 
test  n,  items  until  all  fail. 

•  number  of  component  types,  k  :  5   and  15 

•  system  reliability,  Rs  :  0.9   and  0.975 

•  significance  level,  a.  :  0.05  and  0.2 

•  component  time,  /,  :  small  to  large  (see  tables) 

•  sample  size,  //,  :  small  to  large  (see  tables) 

k 
Reliability  of  a  series  system  is  expressed  as  Rs  =    exp{  —  £  //,}  .   We  chose  arbitrarily 

to  determine  the  failure  rates,  /,  ,  from  this  equation  by  assuming  all  /,r,  to  be  equal, 

consequently. 


-\nRc 

c.     Replications 

The  procedure  was  replicated  1000  times  for  each  case  to  get  1000  values 
°f  Rs,lm  and  1000  values  of  RSMt)  .  We  order  each  set  of  the  1000  values  of  RSM%)  and 
Rs.lm  m  ascending  order.  Then  the  two  approximate  confidence  bounds.  /£s,r.w<ioooa-«»  anc* 
Rs,lm  doood  -»»  -  for  &s  are  tne  1000(1  —  a)th  order  values  of  these  two  sets  of  data.  Finally, 
for  each  of  the  two  ordered  sets  of  data,  we  find  the  order  indices  j\  and  j2  for  which 
Rs,L(*)U\)  an^  ^s.z.(«)(;2)  are  cl°sest  t0  R-s  for  their  respective  sets.  Then  j}  j  1000  and 
j2  j  1000  are  called  the  two  corresponding  simulated  true  confidence  levels. 
2.     Results 

Tables  1  through  3  show  the  results  of  the  3  cases  simulated.  The  results  indi- 
cate that  this  interval  estimation  method  using  estimates  of  failure  rate  ratios  will  yield 
quite  accurate  lower  confidence  limits  for  system  reliability  when  components  have  un- 
known constant  failure  rates.  The  jacknife  method  also  yields  very  accurate  confidence 
limits  which  are  slightly  conservative.  That  is  the  100(1  —  a)  percentile  points  of  RSLM  , 
given  in  the  tables,  are  slightly  less  than  the  true  value  of  Rs.  Consquently, 
P{Rs,L{x)  ^  Rs)  >  1  —  a  .  Alternatively,  one  can  say  RSMa)  is  a  conservative  100(1  —  a) 
percent  lower  confidence  limit  procedure  for  Rs  . 


Table   1.     RELIABILITY   OF  A  SERIES   SYSTEM   WITH   SMALL   NUMBER 
(  LESS  THAN  10  )  OF  SAMPLE  SIZES  -  EXPONENTIAL  CASE 


Number  of 
Compo- 
nents 

Reliability 
of  System 

a 

Lower  Confidence 

Limit 

True  Confidence 
Limit 

WO 

Jacknifing 

WITH 

Jacknifing 

WO 

Jacknifing 

WITH 

Jacknifing 

5 

.90 

2 

0.9017 

0.8904 

0.7720 

0.9 100 

.05 

0.8985 

0.8889 

0.9620 

0.9780 

.975 

2 

0.9750 

0.9720 

0.8020 

0.9160 

.05 

0.9746 

0.9729 

0.9600 

0.9800 

15 

.90 

2 

0.9019 

0.8901 

0.7600 

0.9500 

.05 

0.9012 

0.8933 

0.9300 

0.9800 

.975 

-> 

0.9745 

0.9724 

0.8300 

0.96DO 

.05 

0.9746 

0.9725 

0.9600 

0.9  800 

Sample  sizes  for  5  components  are  9,  7.  10,  8,  6 

Sample  sizes  for  15  components  are  6,  7.  5.  6,  7,  8.  9.  5.  S.  7,  9.  10.  7,  9.  6 

Component  times  /,  are  2.  3.  7.  8.  10,  3,  7,  10.  1,  7.  8,  3.  10,  1,  8 


Table  2.     RELIABILITY  OF  A  SERIES  SYSTEM  WITH  MEDIUM  NUMBER 
(  LESS  THAN  30  )  OF  SAMPLE  SIZES  -  EXPONENTIAL  CASE 


Number  of 
Compo- 
nents 

Reliability 
of  System 

a 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

WO 

Jacknifing 

WITH 

Jacknifing 

W  O 

Jacknifing 

WITH 

Jacknifing 

5 

.90 

2 

0.8999 

0.8978 

0.8020 

0.8660 

.05 

0.8990 

0.89  S3 

0.96(H) 

0.9680 

.975 

2 

0.9750 

0.9745 

0.8000 

0.8460 

.05 

0.9748 

0.9745 

0.9620 

0.9740 

15 

.90 

2 

0.8992 

0.S96O 

0.S3O0 

0.9100 

.05 

0.8980 

0.8981 

0.9  6oo 

0.9700 

.975 

.2 

0.9751 

0.9738 

0.7900 

0.9000 

.05 

0.9750 

0.9741 

0.9500 

0.9800 

Sample  sizes  for  5  components  are  30, 
Sample  sizes  for  15  components  are  3, 


21,  10,  15.  26 

7,  10,  15,  20,  15,  7,  5,  30,  20,  10,  7,  13,  21,  30 


Table  3.     RELIABILITY   OF   A   SERIES   SYSTEM   WITH   LARGE   NUMBER 
(  LESS  THAN  100  )  OF  SAMPLE  SIZES  -  EXPONENTIAL  CASE 


Number  of 
Compo- 
nents 

Reliability 
of  System 

a 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

WO 

Jacknifing 

WITH 

Jacknifing 

WO 

Jacknifing 

WITH 

Jacknifing 

5 

.90 

2 

0.8999 

O.S991 

0.8080 

0.8  2oo 

.05 

0.8990 

0.8992 

0.9680 

0.9620 

.975 

2 

0.9750 

0.9750 

0.8540 

0.8060 

.05 

0.9751 

0.9750 

0.9380 

0.9520 

15 

.90 

2 

0.9000 

0.898S 

0.801O 

0.8750 

.05 

0.9000 

0.8996 

0.9500 

0.9590 

.975 

2 

0.9750 

0.9747 

0.8460 

0.8730 

.05 

0.9750 

0.9750 

0.952" 

0.9470 

Sample  sizes  for  5  components  are  30.  63,  75.  98.  26 

Sample  sizes  for  15  components  are  15.  40.  35.  17,  26.  67,  50.  65.  80.  32.  95.  100,  15, 

45.  30 
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III.     INTERVAL  ESTIMATION  PROCEDURE  -  WEIBULL  CASE 

A.     LOWER  CONFIDENCE  LIMIT   /?5i 

Consider  a  series  system  with  k  components.    Let  the  time  to  failure  X,  of  compo- 
nent /  have  a  IVeibull  distribution  with  density 


ftW   =   Af'fitf'-1  expf-^r/'}     ,        /,>0.  (3.1) 


Then 


Rt{tt)    =     exp{ -(;,•//}  ,         tt>0  (3.2) 


and 


Rs(t)  =  expj  -  V  ;f  /f  1  =  expj  -  SmT  rt  tf<  I  ,       r  >  0  .  (3.3) 

where  /'  =  Af>  ,  /.'„  =  max/'  and  rt  =  ).)\)*m.  If  the  ft,  are  known,  A"f;  has  constant 
failure  rate  /f  and  the  procedures  in  Chapter  II  can  be  used  to  obtain  Rs,Lw  ar*d 
RSMa)  with  T,j  replaced  by  Tfy  in  equation  (2.5).  It  is  well  known  that  }'=lnX,  has 
an  extreme  value  distribution  with  CDF 

Fy{v)    =    l-exp{-<fT~} 

where    t .  =  ln(l//,)   and   <5,=  lift,. 

Engelhardt  and  Bain  [Ref.  17:  p.  323]  have  developed  the  following  simple  unbiased 
estimators  for   c,  and  S,  ,  using  ordered  values     Ytj  =  In  Xm 


di  =  T  =  71 (3-4) 

where  5  =  [0.84/;,]  =  largest  integer  <  0.84  n,   and  Xt(J)  is  the  jlh  order  statistics  from  the 
sample  of  size   ni   of  Xt  .   Also 
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l  =  \n(-L)=yi  +  yd  (3.5) 


where   y  —  0.5772   and  1v,  =  Z^,/^-    Let 


h    =    4"  (3-6) 

and 

Ty    =    Xl'         i  =   1,2,...., 7,      j  =  1,2,...,*.  (3.7) 

We  approximate  the  distribution  of  T:,   by  the  exponential  distribution  with  failure  rate 
)h  =  a'   and  proceed  as  in  Chapter  II.    Define 

*  n-, 

,,     =     -4r  (3.8) 


/:■ 


where     T,    =    Vfi;       /  =   1,2,  •••,/:.      Let   /^  =  max  /'     and 


'■m 


(3.9) 


Then  an  approximate  upper  confidence  limit  for  /'m  is  given  by 


2 
A*  y?.  In 

>-m.im  =   — r —  <3-10) 

£  fa 

(=1 

and  the  corresponding  approximate  lower  confidence  limit   Rs,lm  f°r   ^s(J)   is  given  by 

*SZ(a)     =      eXP  )   "  Vi/(*)X  ^  f  °' ]  l) 

where   n  =  ^n,  . 
i=i 
This  procedure  is  labeled  the  Formula  procedure  in  the  tables  that  follow  in  this 

section.  Its  distinguishing  feature  is  the  equation  for  /?,  given  by  equation  (3.6).  An  al- 
ternative procedure,  labeled  the  Newton  -  Raphson  procedure,  estimates  /?,  using  the 
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maximum  likelihood  procedure  which  is  obtained  using  a  Newton  -  Raphson  approxi- 
mation method.  Equations  for  /?,  under  the  Newton  -  Raphson  procedure  are  provided 
in  Appendix  A. 

B.     SIMULATIONS  AND  RESULTS 
1.     Simulation 

a.  Simulation  language  and  package 

The  programming  language  used  to  simulate  this  problem  was  VS 
FORTRAN  on  an  IBM  3033.  Also  LEXPN  in  LLRANII  was  used  to  generate  observed 
exponential  random  variates.    SHSORT  was  used  to  perform  the  sort  routines. 

b.  Cases  and  Input  parameters 

The  six  input  parameters  below  determine  the  conditions  for  exercising  the 
simulation  runs.  System  reliability  is  determined  by  five  parameters.  Values  of  these 
parameters  are  given  in  the  tables  that  show  the  results.  The  test  plan  simulated  was  to 
test  //,  items  until  all  fail. 

•  number  of  component  types,  k  :  5   and  15 

•  system  reliability,  Rs  :  0.9   and  0.975 

•  significance  level,  a  :  0.05  and  0.2 

•  component  time,  /,  :  small  to  large  (see  tables) 

•  sample  size,  n,  :  small  to  large  (see  tables) 

k 

Reliability  of  a  series  system  is  expressed  as  Rs  =  exp{  —  £  /,/.}  .  Failure  rates  /,  can 
be  determined  from  that  equation  by  assuming   all  /,/,  to  be  equal,  consequently, 

-  In  RK 


A>=  — 


c.     Replications 

The  procedure  was  replicated  1000  times  for  each  case  to  get  1000  values 
of  RSLM  and  1000  values  of  RSMa)  .  We  order  each  set  of  the  1000  values  of  RSiLM  and 
Rs.lw  m  ascending  order.  Then  the  two  approximate  confidence  bounds,  Rs,lmwooo-*»  an^ 
R-s.lm  doood-,))  >  f°r  Rs  are  tne  1000(1  —  a)th  order  values  of  these  two  sets  of  data.  Finally, 
for  each  of  the  two  ordered  sets  of  data,  we  find  the  order  indices  j\  and  j2  for  which 

Rs.lwui)  an<^  Rs,lm(j2)  are  ^e  cl°sest  t0  Rs  f°r  their  respective  sets.  Then  jl  /  1000  and 
j2 1  1000  are  called  the  two  corresponding  true  confidence  levels. 
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2.     Results 

Tables  4  and  5  display  the  results  of  the  simulations  and  determine  the  accuracy  of 
Rs,lw  given  in  equation  (3.11)  as  a  lower  confidence  limit  procedure  for  system  reliability, 
Rs,  for  parameter  values  /„  /?,-, /,  and  Rs  given  in  the  tables.  The  terms  in  the  tables  have 
the  same  meaning  as  the  corresponding  terms  in  Tables  1  through  3  which  were  ex- 
plained in  Section  2.C.2  .  The  procedure  would  be  exact  for  the  Formula  method  if  the 
values  in  the  Formula  column  equal  the  corresponding  numbers  in  the  same  row  in  the 
Reliability  of  System  column.  In  Table  4  for  example,  the  80th  percentile  point  of  .9205 
for  the  Newton  -  Raphson  procedure  is  more  accurate  than  the  Formula  procedure 
which  has  an  80th  percentile  point  of  .9324.  The  last  column  shows,  however,  that  what 
we  have  called  an  80%  lower  confidence  limit  procedure  is  in  fact  closer  to  a  61%  pro- 
cedure. The  accuracy  improves  if  the  sample  size  is  increased  from  15  to  30  as  indicated 
in  Table  5.  The  accuracy  improves  even  more  if  the  number  of  components  in  the  system 
increases  from  5  to  15  as  indicated  in  Tables  4  and  5. 

The  results  in  Tables  4  and  5  show  that  the  RLW  method  given  by  equation 
(3.11)  is  too  optimistic  -  especially  for  small  sample  sizes  and  for  systems  with  a  small 
number  of  components.  Modifications  to  this  procedure  are  needed.  These  tables  also 
indicate  that  the  Newton  -  Raphson  method  is  superior  to  the  Formula  method. 
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Table  A.     ACCURACY  OF  RSL  AS  A  LOWER  CONFIDENCE  LIMIT  FOR  Rs 
WITH  SAMPLE  SIZES  OF  15  -  WEIBl      '   CASE 


Number  of 
Compo- 
nents 

Reliability 
of  System 

a. 

Lower  Confluence 
Limit 

True  Confidence 
Limit 

Formula 

Newton 
Raphson 

Formula 

Newton 
Raphson 

5 

.90 

2 

0.9324 

0.9205 

0.4820 

0.6060 

.05 

0.9512 

0.9406 

0.5380 

0.6660 

.975 

2 

0.9849 

0.9824 

0.5040 

0.6540 

.05 

0.9905 

0.9S91 

0.5600 

0.6660 

15 

.90 

.2 

0.9193 

0.9030 

0.5440 

0.7620 

.05 

0.9416 

0.9288 

0.6060 

0.7760 

.975 

•) 

0.9802 

0.9747 

0.6160 

0.8080 

.05 

0.9S65 

0.9S2S 

0.6580 

0.8240 

Initial  beta  for  5  components  are  2.2.  2.4.  2.6,  2.8.  2.5 

Initial  beta  for  15  components  are  2.2.  2.4,  2.6.  2.8,  2.5,  2.3,  2.7,   2.5,  2.2.  2.6.  2.9,  2.8, 


i  j 


2.6 


1 


Component  times  t,  are  2,  3,  7,  8,  10,  3,  7,  10,  1,  7,  8,  3,  10,  1,  8 


Table  5.     ACCURACY  OF  RSL  AS  A  LOWER  CONFIDENCE  LIMIT  FOR  R5 
WITH  SAMPLE  SIZES  OF  30  -  WEIBULL  CASE 


Number  of 

Compo- 
nents 

Reliability 
of  System 

a. 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

Formula 

Newton 
Raphson 

Formula 

Newton 
Raphson 

5 

.90 

~) 

0.9249 

0.9125 

0.4100 

<>.632o 

.05 

0.9424 

0.0321 

0.5100 

0.6940 

.975 

2 

0.9S37 

0.9798 

O.4380 

0.6160 

.05 

0.9879 

0.9845 

0.4600 

0.6760 

15 

.90 

2 

0.9207 

0.9036 

0.5560 

0.7200 

.05 

0.9333 

0.9198 

0.4  3  SO 

0.7520 

.975 

2 

0.9808 

0.9758 

0.5120 

0.7520 

.05 

0.9850 

0.9810 

0.5200 

0.7420 

Initial  beta  for  5  components  are  2.2,  2.4,  2.6,  2.8,  2.5 

Initial  beta  for  15  components  are  2.2,  2.4,  2.6,  2.8,  2.5,  2.3,  2.7,    2.5,  2.2,  2.6,  2.9,  2.8, 

2.4,  2.6.  2.1 
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IV.     INTERVAL  ESTIMATION  PROCEDURE  -  NORMAL  CASE 

A.     BACKGROUND 

In  recent  years  it  has  become  popular  to  model  mechanical  reliability  as  P(X>  Y)  where 
X  denotes  strength  and  Y  denotes  stress.  This  formulation  of  reliability  is  an  average  of 
P(X>y\  Y  =  y)  ,  since 


P{X>  Y)  =  EyP{X>y)  = 


P{X>y)L(y)dy 


An  alternative  model  is  the  worst  case  approach.  In  the  worst  case  model,  one  se- 
lects a  worst  case  value  ofj\  say  y0  ,  and  then  designs  the  strength.  X  , of  the  component 
so  that  P(X>y0)  =  R0  where  R0  is  a  reliability  requirement.  If  A'  has  a  normal  distrib- 
ution then  this  requirement  imposes  constraints  on  the  mean  and  variance  of  X  .  This 
is  usually  done  in  a  manner  to  comply  with  standard  "safety  factor"  procedures. 

The  average  model,  P(X>  Y),  uses  two  random  variables  and  is  subject  to  more 
random  error  than  the  worst  case  model.  The  accuracy  of  the  expression  P(A'>  }")  for 
values  of  this  expression  close  to  one  is  highly  suspect  due  to  deviations  in  the  tail 
probabilities  of  both  A'  and  }'  from  those  assumed  in  the  model.  It  is  common  to  assume 
that  both  A'  and  }'  have  normal  distributions.  Truncated  normal  distributions  would  be 
more  appropriate  for  many  types  of  mechanical  equipment.  Harris  and  Soms  [Ref.  18: 
pp.  650-663]  discuss  implications  of  this  problem.  Very  significant  errors  in  point  and 
interval  estimation  for  reliability  are  readily  demonstrated  when  X  is  truncated  normal 
but  assumed  to  be  normal  in  the  more  simple  model  which  specifies  that  R  =  P(X  >  }"). 
Table  10  shows  this  effect  when  Zis  truncated  above  at  fj.  +  1.645a  ,  where  Z%  is  the  100 
(1  —  oc)th  percentile  point  of  the  standard  normal  distribution.  Church  and  Harris  [Ref. 
12:  pp.  49-54]  and  Downton  [Ref.  13:  pp.  551-558]  have  developed  approximate  confi- 
dence intervals  for  P(X  >  Y)  when  X  is  normal  with  unknown  mean  and  variance  and  Y 
has  the  standard  normal  distribution. 
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B.     EXACT  AND  APPROXIMATE  INTERVAL  ESTIMATES 

Minimum  variance  unbiased  estimators  (MVUE)  for  R  =  P(X > y)  are  well  known 
when  X  is  normal  with  unknown  mean  and  known  variance  and  also  when  the  variance 
is  unknown.  In  the  former  case,  Lieberman  and  ResnikofT  [Ref.  8]  developed  the  result 


which  is  MVUE  for  R  where  <$>  is  the  standard  normal  cumulative  distribution  function. 
When  the  variance  is  unknown,  several  versions  of  an  integral  expression  for  P(X>y) 
have  been  developed  by  Lieberman  and  ResnikofT  [Ref.  8]  ,  Basu  [Ref.  19]  and  Folks  and 
others  [Ref.  9]. 

Exact  lower  confidence  interval  estimates  for  P{X  >  y)  when  X  is  normal  with  un- 
known mean  and  variance  involve  the  non-central  t  -  distribution.  Owen  and  Ilua  [Ref. 
11]  developed  tables  of  the  lower  90%  and  95%  confidence  limit  values  RL  for  P(X>y) 
based  on  the  non-central  /  -  distribution.  These  values  are  tabled  for  values  of  A'  in  the 
range  -3.0  (.2)  6.0  and  sample  sizes  n  =  2(1)  18.  21  (3)  30.  40  (20)  100  ,  where 
k  =  (x  - y)ls  and  3c  and  5  are  the  sample  mean  and  standard  deviation.  Their  tables  are 
essentially  exact.  An  approximation  to  their  exact  tabled  values  is  given  by  R'L  where 
R'i  =  <D(y')  , 


It, 


(4.1) 


k  =  (x  —y)js  and  fs  „_,  is  the  100(1  —  a)th  percentile  of  the  /  distribution  with  n  —  1  de- 
grees of  freedom.  Equation  (4.1)  was  developed  in  this  thesis.  It  is  an  extensive  ad  hoc 
modification  of  an  equation  developed  by  Church  and  Harris  [Ref.  12].  Tables  6  and  7 
display  values  of  RL  =  O(y)  ,  y  ,  y\  R'L  ,  and  R[  —  RL  for  k  =  1.2,  3.  4.  sample  sizes  n 
=  10.  IS.  30  and  confidence  levels  90%  and  95%  .  RL  and  y  denote  the  "exact"  lower 
confidence  limits  and  corresponding  <t>'l(RL)  values  from  Owen  and  Hua  [Ref.  1 1].  Both 
y'  and  R'L  are  given  by  equation  (4.1).  The  accuracy  of  the  approximate  confidence  in- 
terval is  quite  good  relative  to  the  values  for  RL  given  by  Owen  and  Hua  [Ref.  1 1]. 
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Table  6.     APPROXIMATE  (  R*L  ) 
LIMITS    FOR  P(X>  Y0) 


AND     EXACT  (  RL  )  90%  CONFIDENCE 


n 

*  I0,n-1 

k 

Rvtyy) 

1 

> 

/?[:0(/) 

R]  -  RL 

10 

1.3830 

1 

0.68156 

0.47194 

O.45454 

0.67528 

-0. 00628 

2 

0.8913d 

1.23397 

1.20199 

0.88517 

-0.00613 

0.97453 

1.95262 

1.8S991 

0.97025 

-0.00428 

4 

0.99602 

2.65302 

2.54950 

0.99437 

-0.00165 

18 

1.3334 

1 

0.73037 

0.61385 

0.61133 

0.72950 

-0.00087 

2 

0.92569 

1.44512 

1.44038 

0.92488 

-0.00081 

3 

0.98755 

2.24314 

2.23150 

0.98684 

-0.00071 

4 

0.99877 

3.02679 

3.00614 

0.99859 

-0.00018 

30 

1.3114 

1 

0.75937 

0.70424 

0.70508 

0.75960 

0.00023 

2 

0.94256 

1.57740 

1.57852 

0.94248 

-  0.O000  8 

•3 

0.99233 

2.42404 

2.42459 

0.99205 

-0.00028 

4 

0.99945 

3.2(007 

3.25926 

0.99940 

-  0.0O0O  5 

Table  7.     APPROXIMATE  (   R*L  ) 
LIMITS    FOR  P(X>  Y0) 


AND     EXACT  (  RL  )  95%   CONFIDENCE 


n 

'.OS.n-l 

k 

RL:0(y) 

/ 

R'L:0(f) 

Rl  -  RL 

10 

1.8331 

1 

0.63052 

0.33311 

o.:7-d: 

0.60912 

-0.02140 

■~> 

O.S51S7 

1.04477 

0.94228 

0.82692 

-0.02495 

0.95565 

1.70308 

1.52864 

0.93655 

-0.01910 

4 

0.99031 

2.33813 

2.07743 

0.98075 

-0.00956 

18 

1.7396 

1 

0.69504 

0.51007 

0.49292 

n.(Avx, 

-0.00608 

2 

0.90348 

1.30221 

1.26991 

0.89777 

-0.00571 

■5 

i 

0.97996 

2.05344 

1.99739 

0.97674 

-0. 00322 

4 

0.99735 

2.7S717 

2.70338 

0.99638 

-0.00097 

30 

1.6991 

1 

0.73361 

0.62369 

0.61789 

0.73167 

-0.00194 

2 

0.92855 

1.46579 

1.45391 

0.92677 

-0.00178 

i 
j 

0.98855 

2.27522 

2.25448 

0.9875S 

-0.00097 

4 

0.99894 

3.07140 

3.04028 

0.99873 

-0. 000  21 
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Tables  8  and  9  display  the  results  of  computer  simulations  with  1,000  replications 
to  check  the  accuracy  of  the  R'L  method  for  80%  and  90%  lower  confidence  limits  for 
P{X>y)  s  R  for  j'  =  3,  with  various  values  of  a  and  \x  determined  so  that  R  equals 
the  values  shown.  The  procedure  would  be  exact  at  the  80%  level  if  the  values  in  the 
column  labelled  a  =  .2  equal  the  corresponding  values  of  R  in  the  same  row.  The  "true" 
confidence  level  corresponds  to  the  index  i{R)  of  1,000  ordered  values  of  R'LW  for  which 
R'l(,mr)  =  ^  •  F°r  example,  in  the  seventh  row  of  Table  8,  R  =  .950  ,  ox  =  20  ,  X  =  10, 
R'uwm  was  -9575'  and  R-iawm  was  -9533.  Also  /?l<.20),-5g  =  .950  and  R'Ui0)m  =  .950. 
Tables  8  and  9  indicate  that  the  R'L  procedure  given  by  equation  (4.1)  is  quite  accurate 
at  the  80%  level  of  confidence  for  the  cases  simulated. 
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Table  8.     ACCURACY  ANALYSIS  OF  /?;  PROCEDU 
LOWER     CONFIDENCE     LIMITS     FOR 
NORMALLY  DISTRIBUTED 


RE  FOR  80%  AND  90% 
P(X>3)    WHEN     X     IS 


R 

°x 

n 

'VL.l  000(1-*) 

True  Confidence  Level 

a  =  .2 

a  ---■   .1 

a  =  .2 

a  =  .1 

.950 

0.5 

10 

0.9624 

0.9504 

0.73  lo 

0.8970 

25 

0.9546 

0.9542 

0.7560 

0.8820 

75 

0.9517 

0.9501 

0.7800 

0.89  SO 

1.0 

10 

0.9606 

0.9531 

0.7470 

0.8910 

25 

0.9501 

0.9552 

0.7970 

0.8770 

75 

0.9509 

0.9528 

0.7900 

0.8770 

20.0 

10 

0.9575 

0.9533 

0.7580 

0.8860 

25 

0.9558 

0.9551 

0.74  lo 

0.8730 

75 

0.9511 

0.9511 

0.7830 

0.SS20 

.990 

0.5 

10 

0.9924 

0.9913 

1 .0000 

0.8  840 

25 

0.99<)6 

0.9S98 

0.7800 

0.907( ) 

75 

0.9903 

O.9907 

0.9800 

0.8780 

1.0 

10 

0.9921 

0.9S99 

0.9990 

0.9000 

25 

0.9913 

0.9904 

0.7570 

0.8910 

75 

0.9904 

0.9898 

0.7800 

0.9060 

20.0 

10 

0.9923 

0.9922 

0.9S80 

0.8750 

25 

0.9909 

0.9908 

0.7760 

0.8  Son 

75 

0.9902 

0.9903 

0.7940 

0.8920 

.995 

0.5 

10 

0.9966 

0.9948 

1.0000 

0.9030 

25 

0.9955 

0.9949 

0.7790 

0.9020 

75 

0.9956 

0.9949 

0.7540 

0.9060 

1.0 

10 

0.9966 

0.9952 

i.oooo 

0.8950 

25 

0.9960 

0.9955 

0.8630 

O.SSSO 

75 

0.9949 

0.9947 

0.8040 

0.9150 

20.0 

10 

0.9970 

0.9951 

1.0000 

1.0000 

25 

0.9960 

0.9957 

0.7720 

0.8700 

15 

0.9951 

0.9948 

0.7940 

0.9140 
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Table  9.  ACCURACY  ANALYSIS  OF  R*L  PROCEDURE  FOR  80%  AND  90% 
LOWER  CONFIDENCE  LIMITS  FOR  P(A>30)  WHEN  X  IS 
NORMALLY  DISTRIBUTED 


R 

ox 

n 

"■L,  1000(1  -*) 

True  Confidence  Level 

a  =  .2 

a  ----   .1 

a  -  .2 

a  =  A 

.950 

0.5 

10 

0.9621 

0.9466 

0.75K) 

0.9000 

25 

0.9530 

0.9519 

0.7680 

0.8890 

75 

0.9507 

0.9491 

0.7880 

0.9070 

1.0 

K) 

0.9605 

0.9528 

0.7480 

0.8910 

25 

0.9497 

0.9549 

O.S030 

0.8810 

75 

0.9509 

0.9527 

0.7910 

0.3780 

20.0 

10 

0.9575 

0.9533 

0.7580 

0.8  860 

25 

0.9558 

0.9551 

0.7410 

0.8730 

75 

0.9511 

0.9511 

0.7830 

0.8820 

.990 

0.5 

10 

0.9923 

0.9911 

1.0000 

0.8860 

25 

0.9900 

0.9890 

0.7990 

0.9 140 

75 

0.9901 

0.9904 

0.7980 

0.8860 

1.0 

10 

0.9921 

0.9899 

0.9990 

0.9310 

25 

0.9911 

0.9902 

0.7620 

0.8950 

75 

0.9904 

0.9S98 

0.7800 

0.9070 

20.0 

10 

0.9923 

0.9922 

0.9880 

0.8750 

25 

0.9909 

0.9908 

0.7760 

0.8  8oo 

75 

0.9902 

0.9903 

0.7940 

0.8840 

.995 

0.5 

10 

0.9965 

0.9946 

1.0000 

0.9040 

25 

0.9952 

0.9946 

0.7900 

0.9140 

75 

0.9955 

0.9948 

0.7620 

0.91O0 

1.0 

10 

0.9966 

0.9951 

1 .0000 

0.9280 

25 

0.9959 

0.9953 

0.76 10 

0.8920 

75 

0.9949 

0.9947 

0.8040 

0.9 140 

20.0 

10 

0.9970 

0.9951 

1 .0000 

1.0000 

25 

0.9960 

0.9957 

0.7600 

0.8700 

75 

0.9951 

0.9948 

0.7940 

0.9340 
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Table  10  displays  the  inaccuracies  of  the  R[  lower  confidence  limit  procedure  for 
P(X>  3)  using  equation  (4.1),  which  assumes  X  is  N(n,  a2),  when  in  fact  X  has  the  dis- 
tribution of  a  normal  random  variable  with  mean  \x  and  variance  o-  that  has  been  trun- 
cated at  a  +  1.645(7.   Note  that  y.  +  1.645a  >  >  3,   because  P{X  >  3)  =  .95,  .99  and  .995. 

Examination  of  Table  10  reveals  gross  inaccuracies;  consequently,  even  when  the 
distribution  of  X  is  truncated  far  above  the  value  y,  the  exact  lower  confidence  limit  for 
P(X>y)  can  be  greatly  in  error  when  computed  under  the  assumption  that  X  is  normal. 
This  problem  will  be  compounded  when  one  is  computing  "exact"  confidence  intervals 
for  P{X>  }*)  assuming  both  X  and  Y  are  normal  when  in  fact  one  or  both  may  be  trun- 
cated normal.  This  suggests  that  modeling  mechanical  reliability  as  P{X>  Y)  may  be 
more  risky  than  more  standard  engineering  approachs  for  modeling  mechanical  reliabil- 
ity which  include  the  notion  of  safetv  factors. 
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Table  10.  ACCURACY  ANALYSIS  OF  R*L  PROCEDURE  FOR  80%  AND  90% 
LOWER  CONFIDENCE  LIMITS  FOR  P(X>3)  WHEN  X  IS  TRUN- 
CATED NORMAL 


Reli- 
abil- 
ity of 
Sys- 
tem 

o, 

n 

'V. 1000(1 -i) 

True  Confidence  Level 

a  =  .2 

a  =  .1 

a.  =  .2 

a  =  .1 

.950 

0.5 

10 

0.9806 

0.0506 

1.0000 

1.0000 

25 

0.8941 

0.0125 

1 .0000 

1 .0000 

75 

0.9383 

0.0871 

1.0000 

1 .00( )( ) 

1.0 

10 

0.S904 

0.0419 

1.0000 

1.0000 

25 

0.9555 

0.0813 

1 .0000 

1.0000 

75 

0.9431 

0.0782 

1.0000 

1.0000 

20.0 

10 

0.8503 

0.0105 

1 .0000 

1.0000 

25 

0.9274 

0.04S5 

1 .0000 

1.0000 

75 

0.9496 

0.1222 

1 .0000 

1 .0000 

.990 

0.5 

10 

0.9734 

0.6035 

1 .0000 

1.0000 

25 

0.9924 

0.6346 

1  .oooo 

1  .oooo 

75 

0.9780 

0.6552 

1.0000 

1.0000 

1.0 

1<> 

0.9976 

0.9172 

1.0000 

1.0000 

25 

0.9S84 

0.7505 

1  .oooo 

I.  oooo 

75 

0.9951 

0.6800 

1.0000 

1 .0000 

20.0 

10 

0.9593 

0.5116 

1.0000 

1 .0000 

25 

0.9989 

0.6821 

1.0000 

1.0000 

15 

0.9SS3 

0.7514 

1.0000 

1 .0000 

.995 

0.5 

10 

0.9813 

0.7328 

1  .Oooo 

1.0000 

25 

0.9719 

0.9096 

1.0000 

1.0000 

75 

0.9933 

0.9130 

l  .oooo 

1.0000 

1.0 

10 

0.9709 

0.9015 

1.0000 

1.0000 

25 

0.9991 

0.8996 

1.0000 

1.0000 

75 

0.9949 

0.9267 

1.0000 

1 .0000 

20.0 

10 

0.9723 

0.6432 

1.0000 

1.0000 

25 

0.9923 

0.7539 

1.0000 

1 .0000 

75 

0.9971 

0.9489 

1.0000 

1.0000 
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V.     CONCLUSIONS  AND  RECOMMENDATIONS 

The  lower  interval  estimation  procedure  for  reliability  of  coherent  systems  which 
was  developed  in  Chapter  II  appears  to  be  accurate,  easy  to  use,  and  applicable  to  co- 
herent systems.  Although  this  procedure  assumes  that  failure  times  of  all  components 
of  the  system  are  independent  and  have  exponential  probability  distributions,  it  can  be 
easily  extended  to  systems  with  component  failure  distributions  that  can  be  transformed 
into  exponential  failure  distributions  ;  e.g.,  U'eibull  distribution  with  known  shape  pa- 
rameter. This  procedure  has  potential  for  being  combined  with  a  similar  procedure  for 
systems  with  cyclical  components.  The  combined  procedure  would  provide  for  use  of 
binomial  component  test  data  and  exponential  component  test  data  to  compute  lower 
confidence  limits  on  the  reliability  of  coherent  systems  with  both  cyclic  and  continuous 
components.  Such  an  extension  could  use  a  failure  rate  ratio  estimation  procedure  sim- 
ilar to  that  developed  in  this  thesis.  Such  a  method  should  be  explored. 

The  interval  estimation  method  for  the  reliability  of  a  system  with  components  that 
have  U'eibull  failure  times  is  not  sufficiently  accurate  to  be  applied  to  systems  that  have 
10  or  fewer  components  each  with  ten  or  fewer  tests.  This  procedure  needs  further  study 
and  refinement. 

The  approximate  lower  confidence  limit  for  component  reliability  P(X>y)  when  A' 
is  normally  distributed  with  unknown  mean  and  variance  is  very  accurate.  It  does  not 
require  an  extensive  set  of  tables  such  as  those  developed  by  Owen  and  HuafRef.  11). 
but  only  requires  the  use  of  the  standard  /  tables. 
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APPENDIX  A.     MLE  OF  VVEIBULL  PARAMETERS  BY  THE  NEWTON  - 

RAPHSON  METHOD 

Let  X  ~  WEI(cl,  /?).  Then  the  likelihood  function  for  the  first  r  ordered  observations 
from  a  random  sample  of  size  n  is  given  by 


f{xx  .„,•••,  xr,  n)  = 


(«  "  r)\ 


n\ 


{n  -  r)\ 


i    \    y. 


\fxWn) 

"fl 


u  -  F,{xr:„)y 


-v_1  «p{-[Z(ila-)',+(-<^L)' 


R 
J 


Setting  the  partial  derivatives  of  this  likelihood  with  respect  to  a  and  /?  equal  to  zero 
gives  the  MLE's  a  and  /?  as  solutions  to  the  equations 


and 


/    xf  In  .x 


n 

1  y 


n 


In  X: 


i=\ 


i=\ 


r  V  J 


a.  = 


i=\ 


where  n  =  r.  It  can  be  shown  that  these  equations  have  unique  solutions  which  are  the 
maximum  likelihood  estimates.  The  NEWTON  -  RAPHSON  method  for  solving  an 
equation  g(j$)  =  0  is  to  determine  successive  approximations  /?,,  where 
/L,  =  /),  -  g(fi:)!g'{P,)-   Therefore,  the  estimates  of  a  and  /?  can  be  solved  by  letting 
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gm  = 


/=] 


xf  In  xt 


Y  J 


1  Vi 


The  derivative  of  #(/?)  is 


2>f( 


In  *,-) 


YJ 


I  ln-Y/ 


*'(/*)  = 


/=1 


/=1 


+ 


YJ 


YJ 


i=\ 


i=\ 
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APPENDIX  B.     FORTRAN  CODE  FOR  INTERVAL  ESTIMATION 
PROCEDURE  -  EXPONENTIAL  CASE 

PROGRAM  EXPONE 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAiVAAAAAAAAAAAAA 


c 

c 

* 

c 

A 

c 

-A- 

c 

A 

c 

A 

c 

A 

c 

Vc 

c 

A 

c 

A 

c 

* 

c 

c 

* 

c 

Vc 

c 

-.V 

c 

•>v 

c 

-V 

c 

A 

c 

* 

c 

A 

c 

* 

c 

•J- 

c 

A 

c 

~v 

c 

* 

c 

* 

c 

* 

c 

.'. 

c 

* 

c 

J. 

c 

* 

c 

A 

c 

A 

c 

,•. 

c 

* 

c 

* 

c 

* 

c 

A 

c 

A 

c 

A 

c 

A 

c 

* 

c 

* 

c 

i% 

c 

A 

c 

Jm 

c 

A 

c 

J" 

THIS  IS  A  PROGRAM  TO  COMPUTE  THE  TRUE  CONFIDENCE  LIMIT  OF 
A  SERIES  SYSTEM  WITH  /  WITHOUT  JACKNIFE  METHOD  AND 
COMPARE  THE  DIFFERENCE  OF  THOSE  RESULTS. 

BELOWS  ARE  GIVEN  OR  ASSUMED. 


N 

RSYS 
ALPHA 
TIME 


NUMBER  OF  COMPONENTS  TYPE 
RELIABILITY  OF  A  SERIES  SYSTEM 
SIGNIFICANCE  LEVEL 
TEST  TIME 


SAMPLE  ;  SAMPLE  SIZE  FOR  EACH  TYPE  OF  COMPONENT 


THESE  ARE  VARIABLES  USED. 


B 

BIGLAM 

CASE 

CHISQR 

CLR1 

CLR2 

CTIME 

DF 

DIFF 

LAMBDA 

LAMHAT 

LAMHM 

LAMHST 

LAMMAX 

LOMIT 

KEY1 

KEY2 

R 

Rl 

R2 

RATIOl 

RATI02 

ROM  IT 

ROMSUM 

RSTAR 

RVAL1 

RVAL2 

T 

ZALPHA 


TEMPORARY  ARRAY  FOR  EXPONENTIAL  RANDOM  VARIATE 

LARGEST  VALUE  OF  LAMBDA 

COUNTER 

CHI -SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF 

TRUE  CONFIDENCE  LIMIT  FROM  NON- JACKNIFING 

TRUE  CONFIDENCE  LIMIT  FROM  JACKNIFING 

TEST  TIME  OF  INDIVIDUAL  COMPONENT 

DEGREE  OF  FREEDOM  FOR  CHI -SQUARE. 

ABSOLUTE  VALUE  OF  DIFFERENCE  BETWEEN  Rl  &  R2 

FAILURE  RATE  OF  EACH  COMPONENT  TYPE 

LAMBDA  HAT  FOR  JACKNIFING 

FINAL  LAMBDA  HAT  FROM  NON- JACKNIFING 

FINAL  LAMBDA  HAT  FROM  JACKNIFING 

LARGEST  VALUE  OF  LAMHAT 

LAMBDA  WITH  1  COMPONENT  OMITTED 

TEMPORARY  ARRAY 

TEMPORARY  ARRAY 

INITIALLY  COMPUTED  RELIABILITY  FOR  EACH 

COMPONENT  TYPE 

COMPUTED  RELIABILITY  BY  NON-JACKNIFING 

COMPUTED  RELIABILITY  BY  JACKNIFING 

RATIO  OF  LAMBDA  FOR  NON-JACKNIFING 

RATIO  OF  LAMBDA  FOR  JACKNIFING 

RATIO  WITH  1  COMPONENT  OMITTED 

USED  FOR  JACKNIFE,  SUM  OF  R  WITH  OMIT  1 

FINALLY  COMPUTED  R  VALUE  (  R  STAR  ) 

R(  500  *  (1-ALPHA)  )  FOR  NON-JACKNIFING 

R(  500  *  (1-ALPHA)  )  FOR  JACKNIFING 

TOTAL  TEST  TIME  OF  EACH  COMPONENT  TYPE 

RIGHT  PERCENTILE  POINT( NORMAL  DISTRIBUTION) 


A 

* 

A 
* 
* 
* 
A 
A 
A 
A 
A 
A 
A 
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c  *  * 

Q  lY^WoWcVrVriViiWeVnVycVrVoWf  ******  AVr?^ 

c 
c 

PARAMETER  (NN  =  500) 

REAL  TIME(15),  ALPHA(2),  ZALPHA(2),  RSYS(2),  LAMBDA(15),  P 

REAL  LAMMAX,  LAMHAT(15),  B(100),  RSTAR(15),  LAMHM 

REAL  RATI01(15),  RATI02(15),  BIGLAM,  DIFF(15) 

REAL  LOMIT(15,100),  CTIME( 15 , 100) ,  ROMIT( 15 , 100) ,  T(15) 

REAL  ROMSUM  ,  SUM1,  SUM2 ,  LAMHST,  Rl(NN),  R2(NN),  KEYl(NN) 

REAL  RVALK15),  RVAL2(15),  CLR1(15),  CLR2(15),  KEY2(NN) 

REAL  DIFFR1(15),  DIFFR2(15),  R(15) 

INTEGER   SAMPLE (15),  N(2),  DF,  CASE 

DATA  N  /  5,0   / 

DATA  RSYS  /  .9,  .9  75  / 

DATA  TIME  /  2,3,7,8,10,3,7,10,1,7,8,3,10,1,8  / 

DATA  SAMPLE  /  30,63,75,98,26,15,7,5,30,20,10,7,13,21,30  / 

DATA  ALPHA  /  . 2 ,  .05  / 

DATA  ISEED  /  1736  / 


C     /*   SUBROUTINE  ZTABLE  COMPUTES  THE  RIGHT  PERCENT  POINT  ZALPHA   */ 
C     /•'•-  FROM  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  */ 

DO  10  I  =  1,  2 

CALL  ZTABLE (ALPHA( I),  ZALPHA(I)) 
10  CONTINUE 


CASE  =  0 

/•'•'  II  IS  INDEX  FOR  N  */ 

DO  150  II  =  1,  2 

/*  COMPUTE  THE  DEGREE  OF  FREEDOM  FOR  CHI -SQUARE  */ 

DF  =  0. 

DO  20  I  =  1,  N(II) 

DF  =  DF  +  SAMPLE(I) 
20      CONTINUE 

DF  =  2  *  DF 


C  /*  FINDING  COMPONENT  TYPE  THAT  HAS  THE  MINIMUM  NUMBER      */ 

C  /*  OF  SAMPLE  SIZE  */ 

NSTAR  =999 

DO  25  I  =  1,  N(II) 

IF  (SAMPLE(I)  .  LE.  NSTAR)  THEN 

NSTAR  =  SAMPLE(I) 
END  IF 
25      CONTINUE 
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/*  J J  IS  INDEX  FOR  RSYS  */ 

DO  140  JJ  =  1,  2 

C  /■•  COMPUTE  LAMBDA  FROM  THE  GIVEN  EQUATION  AND        */ 

C  /*  FIND  THE  BIGGEST  LAMBDA  */ 

BIGLAM  =  0 

DO  30  K  =  1,  N(II) 

LAMBDA(K)  =  (  -ALOG(RSYS( JJ) )  /  N(II))  /  TIME(K) 
IF(  LAMBDA(K)  . GE.  BIGLAM)   BIGLAM  =  LAMBDA(K) 
R(K)  =  EXP(  -  LAMBDA(K)  *  TIME(K)  ) 
30  CONTINUE 


/*  COMPUTE  CHI -SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF    */ 


/*  LL  IS  INDEX  FOR  ALPHA  */ 

DO  130  LL  =  1,  2 

IF  (  DF  . EQ.  1  )  THEN 
P  =  ALPHA(LL)  /  2 
CALL  ZTABLE(P,  ZALPHA(LL)) 
CHISQR  =  ZALPHA(LL)  **  2 

ELSE  IF  (  DF  .EQ.  2  )  THEN 

CHISQR  =  -2  *  ALOG(ALPHA(LL)) 

ELSE  IF  (  DF  .  GE.  3  )  THEN 

CHISQR  =  DF  *  (1  -  2./(9  *  DF)  + 

ZALPHA(LL)  *  SQRT(  2./(9  *  DF)))  **  3 
END  IF 


CASE  =  CASE  +  1 
DIFF(CASE)  =  0 

DO  120  L  =  1,  NN 

LAMMAX  =  -99. 

DO  50  I  =  1,  N(II) 


C  /*  GENERATE  EXPONENTIAL  RANDOM         */ 

C  /*  NUMBERS  WITH  MU  =  1  */ 

CALL  LEXPN(ISEED,  B,  SAMPLE(I),  1,  0) 
T(I)  =  0 

DO  40  J  =  1,  SAMPLE(I) 

C  /*  CONVERT  TO  EXPONENTIAL  RANDOM  */ 

C  /*  NUMBERS  WITH  MU  =  LAMBDA  AND   */ 


2S 


C  /*  ADD  THOSE  FOR  EACH  COMPONENT   */ 

C  /*  TYPE  */ 

B(J)  =  B(J)  /  LAMBDA(I) 
T(I)  =  T(I)  +  B(J) 
CTIME(I,J)  =  B(J) 
40  CONTINUE 

LAMHAT(I)  =  SAMPLE(I)  /  T(I) 


C  /*  FINDING  MAXIMUM  LAMBDA  HAT  AND  ITS  */ 

C  /*  INDEX  */ 

IF  (  LAMHAT(I)  . GE.  LAMMAX  )  THEN 

M  =  I 

LAMMAX  =  LAMHAT(I) 
ENDIF 

50  CONTINUE 


C  /*  RATIOl  IS  FOR  WITHOUT  JACKNIFE  */ 

C  /*  RATI02  IS  FOR  WITH    JACKNIFE  */ 

DO  60  I  =  1,  N(II) 

RATIOl(I)  =  LAMBDA(I)  /  BIGLAM 
RATI02(I)  =  LAMHAT(I)  /  LAMMAX 
60  CONTINUE 


C  /*  PART  OF  JACKNIFE  METHOD  FOR  LAMBDA       */ 

C  /*  WITH  OMIT  1  VARIABLE  EACH  TIME  */ 

DO  90  I  =  1,  N(II) 

DO  80  J  =  1,  SAMPLE(I) 

LOMIT(I,J)  =  0. 

DO  70  K  =  1,  SAMPLE(I) 

IF  (J  .NE.  K)  THEN 

LOMIT(I.J)  =  LOMIT(I,J)  +  CTIME(I,K) 
ENDIF 

70  CONTINUE 

LOMIT(I,J)  =  (SAMPLE(I)-l)  /  LOMIT(I,J) 

80  CONTINUE 

90  CONTINUE 
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C  /*  ADAPT  ABOVE  RESULT  TO  OUR  EQUATION  TO    */ 

C  /*  GET  THE  RELIABILITY  AND  TRUE  CONFIDENCE  */ 

C  /*  LIMIT  */ 

SUM1  =  0. 

SUM2  =  0. 

SUM3  =  0. 

SUM4  =  0. 

DO  110  I  =  1,  N(II) 

C  /*  NON  JACKNIFING  (  ORIGINAL  )  */ 

SUM1  =  SUM1  +  RATIOl(I)  *  T(I) 
SUM2  =  SUM2  +  RATIOl(I)  *  TIME(I) 

C  /*  WITH  JACKNIFING  */ 

ROMSUM  =  0. 

DO  100  J  =  1,  NSTAR 

ROMIT(I,J)  =  LOMIT(I.J)  /  LOMIT(M,J) 
ROMSUM     =  ROMSUM  +  ROMIT(I,J) 
100  CONTINUE 

RSTAR(I)  =  NSTAR  *  RATI02(I)  - 
*  (NSTAR  -  1)  *  ROMSUM  /  NSTAR 

SUM3  =  SUM3  +  RSTAR(I)  *  T(I) 
SUM4  =  SUM4  +  RSTAR(I)  *  TIME(I) 
110  CONTINUE 


C  /*  Rl  ;  RELIABILITY  OF  A  SYSTEM  WITHOUT     */ 

C  /*  JACKNIFE  */ 

C  /*  R2  ;  RELIABILITY  OF  A  SYSTEM  WITH        */ 

C  /*  JACKNIFE  */ 

LAMHM  =  CHISQR  /  (2  *  SUM1) 
R1(L)  =  EXP(  -  LAMHM  *  SUM2  ) 

LAMHST  =  CHISQR  /  (2  *  SUM3) 
R2(L)    =  EXP(  -  LAMHST  *  SUM4  ) 


IF  (  ABS(R1(L)  -  R2(L))  . GE.  DIFF(CASE)  )  THEN 
DIFF(CASE)  =  ABS(R1(L)  -  R2(L)) 
DIFFRl(CASE)  =  R1(L) 
DIFFR2(CASE)  =  R2(L) 

ENDIF 

120  CONTINUE 

C  /*  NONIMSL  LIBRARY  'SHSORT'  WILL  SORT  Rl,  R2     */ 

C  /*  BY  SHELL  SORT  ALGORITHM  */ 
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CALL  SHS0RT(R1,  KEY1,  NN) 
CALL  SHS0RT(R2,  KEY2 ,  NN) 

MM  =  NN  *  (1  -  ALPHA(LL)) 

RVALl(CASE)  =  Rl(MM) 
RVAL2(CASE)  =  R2(MM) 

C  /*  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  Rl,  R2    */ 

C  /*  WHICH  THE  VALUE  OF  IT  IS  CLOSEST  TO  RSYS.      */ 

CALL  FINDJ(R1,  NN,  RSYS(JJ),  Jl) 
CALL  FINDJ(R2,  NN,  RSYS(JJ),  J2) 

CLRl(CASE)  =  Jl  /  FLOAT(NN) 
CLR2(CASE)  =  J2  /  FLOAT(NN) 

130  CONTINUE 

140      CONTINUE 
150  CONTINUE 


WRITE(6,600) 
WRITE(6,650) 
WRITE(6,670) 

CASE  =  1 

DO  210  II  =  1,  2 
DO  200  JJ  =  1,  2 
DO  190  LL  =  1,  2 

WRITE(6,700)  N(II),  RSYS(JJ),  ALPHA(LL),  RVALl(CASE), 
*  RVAL2(CASE),  DIFF(CASE),  CLRl(CASE),  CLR2(CASE) 

WRITE(6,888)  DIFFRl(CASE) ,  DIFFR2(CASE) 
CASE  =  CASE  +  1 

190       CONTINUE 
200    CONTINUE 

WRITE(6,777)  (  SAMPLE(J),  J=1,N(II)  ) 

WRITE(6,999)  (  R(J),  J=1,N(II)  ) 
210  CONTINUE 


600  FORMAT(lll,5(/)|7X,'****  RELIABILITY  OF  SERIES  SYSTEM  ****') 

650  FORMAT(///,T50,fRl,  CL1  ;  WITHOUT  JACKNIFING' , 
/,T50.'R2,  CL2  ;  WITH    JACKNIFING') 

670  FORMAT (///,T6, 'NUMBER  OF1 ,T19 ,' RELIABILITY*  ,T33, 'ALPHA' ,T46, 'Rl' , 
T56/R2' ,T63,'MAX  (Rl  -  R2) ' ,T79 ,' TRUE ' ,T87  ' TRUE ',/ , 
T6, 'COMPONENTS' ,T19,' OF  SYSTEM' ,T79,'C.L.  l',T87,'C.L.  2', 
/,T5,89('-')) 
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700   FORMAT(/,T8,I5,T22,F5.  3,T33,F5. 2 ,T43 ,F7. 4 ,T53 ,F7. 4 ,T65 ,F8.  4 ,T78 , 
*  F7.  4,T86,F7.  4) 

777      FORMAT(/,T3,' SAMPLE   SIZE   ARE    ',    5(2X,I3)    ) 
888      FORMAT(T60,fRl=' ,F6.4,'      R2=' ,F6. 4) 
999      F0RMAT(/,T3,'R(I)    ;     ' ,5(  IX, F5.  3) ) 

STOP 

END 


JLJLJL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  J-  J -  JL  J-  *L  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  -'-  -'-  JL  JL  -'-  -'-  JL  - '  -  -'-  JL  JL  J -  JL  - L  JL  .-'-  JL  JL  JL  JL  JL  -' -  y-  -  L  -. '  -  JL  JL  JL  J-  JL  JL  .V  JL  JL  JL  JL  JL  JL  JL  JL  J-  JL  -L 

*  * 

SUBROUTINE   ZTABLE( ALPHA, Z ALPHA) 
C 

C  «   SUBROUTINE    ZTABL2   COMPUTES   RIGHT   PERCENT   POINT   ZALPHA  FROM      » 

C  «   RIGHT  CUMULATIVE   PROBABILITY  ALPHA  » 

C 

REAL  ALPHA,    ZALPHA 

IF    ((ALPHA   .GT.     0.0)    .OR.     (ALPHA   .  LT.     1))    THEN 

V  =    -   AL0G(4   *   ALPHA  *   (    1    -   ALPHA    )) 

ZALPHA  =   SQRT(W  *    (2.06118    -    (    5.72622    /    (W  +    11.6406)))) 

IF    (ALPHA    .GT.     0.  5)    THEN 

ZALPHA  =   -    ZALPHA 
END  IF 
END  IF 

RETURN 
END 

JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  -J-  JL  J.  JL  JL  JL  JL  JL  JL  »L  .J.  JL  JL  JL  JL  J,  .y  J^  JL  J-  J-  -f-  J*  JL  J.  JL  J.  JL  J-  J.  J*  J.  J.  JL  J.  J*  JL  JL  JL  J*  J-  J/,  JL  JL  JL  .L  J.  JL  JL  JL  JL  JL  JL  jl  JL  JL  j.  JL  jl  JL  JL 

SUBROUTINE    SHSORT( A,KEY,N) 
DIMENSION  A(N),KEY(N) 
Ml=l 
6   M1=M1*2 

IF    (Ml    . LE.    N)      GO  TO   6 
Ml=Ml/2-l 
MM=MAX0(M1/2,1) 
GO  TO   21 

20  MM=MM/2 

IF    (MM   .LE.     0)      GO  TO    100 

21  K=N-MM 
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22  DO  1  J=1,K 

II=J 
11  IM=II+MM 

IF  (A(IM)  . GE.  A(II))   GO  TO  1 
TEMP=A(II) 
IT=KEY(II) 
A(II)=A(IM) 
KEY(II)=KEY(IM) 
A(IM)=TEMP 
KEY(IM)=IT 
I 1=1 I -MM 

IF  (II  . GT.  0)   GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 


V?  Vr  VoWw  "  V  ?V  Vw  c  Vr  Vr  Vr  iWr  VoV  ?V  }V  iV  iV  ?V  iV  Vc  VoV  ^^ 

SUBROUTINE  FINDJ(A,  NN,  R,  J) 
C 

C     «  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  ARRAY  A  WHICH  THE      » 
C     «  VALUE  OF  IT  IS  CLOSEST  TO  R.  » 


C 


REAL  A(NN),  R,  VALUE 
INTEGER  J 

VALUE  =  ABS(A(NN)  -  R) 

DO  100  I  =  NN-1,  1,  -1 

IF  (ABS(A(I)  -  R)  .  LT.  VALUE)  THEN 
VALUE  =  ABS(A(I)  -  R) 
ELSE 

J  =  I  +  1 
RETURN 

END  IF 

100  CONTINUE 

RETURN 
END 
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APPENDIX  C.     FORTRAN  CODE  FOR  INTERVAL  ESTIMATION 


PROCEDURE  -  WEIBULL  CASE 


PROGRAM  WEIBUL 


c 

c 

ft 

c 

* 

c 

* 

c 

ft 

c 

-V 

c 

a 

c 

* 

c 

-V 

c 

ft 

c 

* 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

* 

c 

ft 

c 

ft 

c 

-.'.- 

c 

ft 

c 

-'.- 

c 

* 

c 

* 

c 

ft 

c 

-'- 

c 

ft 

c 

* 

c 

-V 

c 

ft 

c 

ft 

c 

ft 

c 

* 

c 

* 

c 

* 

c 

ft 

c 

ft 

c 

v.* 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

ft 

c 

•J. 

ftftftftftftftftftftftftftftftftiVftftftftftftftftftftftftftftftftftftftVfftiVftftftiVftftVcVfftftft^VftftftftVfftftVfftftftftVrftft 


THIS  PROGRAM  COMPUTES  THE  RELIABILITY  OF  A  SERIES  SYSTEM 
USING  TWO  DIFFERENT  METHODS  FOR  THE  CASE  WHEN  FAILURE 
OF   EACH  COMPONENTS  IS  WEIBULLY  DISTRIBUTED. 
AND  ALSO  COMPARES  THE  ESTIMATE  OF  SHAPE  PARAMETERS  (  BETA 
WHICH  IS  COMPUTED  FROM  TWO  DIFFERENT  METHODS. 

THE  SERIES  SYSTEM  CONSIDERED  IN  THIS  PORGRAM  HAS 

N  COMPONENT  TYPES  AND  EACH  COMPONENT  HAS  SAME  NUMBER  OF 

SAMPLE  SIZES. 

BELOWS  ARE  GIVEN  OR  ASSUMED  ; 


BETA 
ALPHA 
X 
RSYS 


SHAPE  PARAMETER  OF  WEIBULL  DISTRIBUTION 

SIGNIFICANCE  LEVEL 

TEST  TIME 

RELIABILITY  OF  A  SERIES  SYSTEM 


THESE  ARE  VARIABLES  USED 


BETHAT 

BETNEW 

CHISQR 

CLR 

CLRNEW 

EXPO 

KSUBM 

LAMBDA 

LAMB I G 

LAMHAT 

LAMMAX 

LAMMU 

LAMSTR 

LAMNEW 

LAMUNR 

RATIO 

RATNEW 

RHAT 

RHTNEW 

XBENEW 

XHAT 

XIJHAT 

XIJNEW 

XNEW 

XTOBET 

W,WEIB 


BETA  HAT  FROM  THE  FORMULA 

BETA  HAT  FROM  NEWTON  -  RAPHSON  METHOD 

CHI -SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF 

TRUE  CONFIDENCE  LIMIT  FROM  THE  FORMULA 

TRUE  CONFIDENCE  LIMIT  FROM  THE  N-R  METHOD 

EXPONENTIAL  RANDOM  NUMBERS  WITH  MU=1 

K  SUB  M  VALUE 

SCALE  PARAMETER  OF  WEIBULL  DISTRIBUTION 

MAXIMUM  VALUE  OF  LAMBDA  FROM  THE  N-R  METHOD 

ESTIMATES  OF  LAMBDA  FROM  THE  FORMULA 

MAXIMUM  VALUE  OF  LAMBDA  FROM  THE  FORMULA 

FINALLY  COMPUTED  LAMBDA  FROM  THE  FORMULA 

LAMBDA  TO  THE  BETA 

ESTIMATES  OF  LAMBDA  FROM  THE  N-R  METHOD 

FINALLY  COMPUTED  LAMBDA  FROM  THE  N-R  METHOD 

RATIO  OF  LAMBDA  FOR  THE  FORMULA 

RATIO  OF  LAMBDA  FOR  THE  N-R  METHOD 

R(  500  *  (1 -ALPHA)  )  FOR  THE  FORMULA 

FOR  THE  N-R  METHOD 


) 


R(  500  *  (1 -ALPHA) 
X  TO  THE  BETNEW 
SUM  OF  XIJHAT 
W  TO  THE  BETHAT 
W  TO  THE  BETNEW 
SUM  OF  XIJNEW 
X  TO  THE  BETHAT 
WEIBULL  RANDOM  NUMBERS 


ft 

ft 
ft 
ft 
ft 
ft 


ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 


ft 

ft 

ft 

* 
ft 

ft 


ft 
ft 
ft 
ft 
ft 
ft 
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C     *       ZALPHA  ;  RIGHT  PERCENTILE  POINT  (NORMAL  DISTRIBUTION)     * 
C 

C 

c 

C     PARAMETER  (  N=5 ,  M=30  ,  NN=  500  ) 
PARAMETER  (  N=15,  M=30  ,  NN=  500) 

REAL  BETA(N),  LAMBDA(N),  X(N),  XHAT(N),  XTOBET(N),  XX(N,M) 

REAL  W(N,M),  Y(N,M),  BETHAT(N),  WBHAT(N,M),  LAMHAT(N) 

REAL  RATIO(N),  LAMSTR(N),  LAMMAX,  RSYS(2),  KSUBM(N) ,  LAMMU 

REAL  EXPO(M),  ALPHA(2),  ZALPHA(2),  RHAT(IO),  RSL(NN),  CHISQR 

REAL  B(M),  KEY(M),  KEYl(NN) ,  KEY2(NN),  XIJHAT(N,M) 

REAL  WEIB(M),  BETNEW(N),  XBENEW(N),  XNEW(N),  XIJNEW(N.M) 

REAL  LAMNEW(N),  LAMBIG,  RATNEW(N),  LAMUNR,  RSLNEW(NN) 

REAL  RHTNEW(IO),  CLR(IO),  CLRNEW(IO) 

REAL  BHTBAR,  BNRBAR,  BHTMSE ,  BNRMSE 

INTEGER   S(N),  DF,  I COUNT,  SAMPLE(N),  CASE 

DATA  ISEED   /  1736  / 
DATA  ALPHA   /  .2,  .05  / 
DATA  RSYS    /  .90,  .975  / 

C     DATA  BETA    /  1.2,1.4,1.6,1.8,1.5  / 

C     DATA  BETA  /  1.  2  , 1.  4  , 1.  6  , 1.  8  , 1.  5  , 1.  3  , 1.  7  , 1.  5  , 1.  2  , 1.  6  , 1.  9  , 1.  8  , 1.  4  , 

C    *  1.  6,1.  1/ 

C     DATA  BETA    /  2.2,2.4,2.6,2.8,2.5  / 

DATA  BETA  /2.  2  ,  2.  4  ,2.  6  ,  2.  8  ,2.  5  ,2.  3  ,  2.  7  ,  2.  5  ,2.  2  ,2.  6  ,2.  9  ,2.  8  ,  2.  4, 
2.6,2.1/ 

C     DATA  X      /  2,3,7,8,10/ 

DATA  X  /2, 3, 7, 8, 10, 3, 7, 10, 1,7, 8, 3, 10, 1,8/ 

C     DATA  SAMPLE  /9, 7, 10, 8, 6/ 

C     DATA  SAMPLE  /30 , 21 , 10  ,  15 , 26/ 

C     DATA  SAMPLE  /  30,63,75,98,26  / 

C     DATA  SAMPLE  /6, 7, 5, 6, 7, 8, 9, 5, 8, 7, 9, 10, 7, 9,  6/ 

C     DATA  SAMPLE  /3 , 7 , 10 , 15 , 20 , 15 , 7 ,5 , 30 , 20 , 10 , 7 , 13 , 21 , 30/ 

C     DATA  SAMPLE  / 15 ,40 , 35 , 17 , 26 , 67 ,50 , 65 , 80 , 32 , 95 , 100 , 15 ,45 , 30/ 

DATA  SAMPLE  /  15  *  15  / 

CASE  =  1 

DO  10  I  =  1,  N 

IF  (  SAMPLE(I)  .LE.  15)   KSUBM(I)  =  1.40 

IF  (  SAMPLE(I)  .LE.  30)   KSUBM(I)  =  1.50 
10  CONTINUE 


DO  220  JJ  =  1,  2 

/*  FINDING  THE  LAMBDA  AND  LAMBDA  STAR  FROM  THE  GIVEN  DATA    */ 

DO  20  I  =  1,  N 
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LAMBDA(I)  =  ((  -  ALOG(RSYS(JJ))  /  N)  **  ( 1. /BETA( I ) ) )/  X(I) 
LAMSTR(I)  =  LAMBDA(I)  **  BETA(I) 
S(I)      =0.84  *  SAMPLE(I)  -  0.5 
20   CONTINUE 

C       /••-  SUBROUTINE  ZTABLE  COMPUTES  THE  RIGHT  PERCENT  POINT  ZALPHA  */ 
C       /*  FROM  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  */ 

DO  40  I  =  1,  2 

CALL  ZTABLE(ALPHA(I),  ZALPHA(I)) 
40   CONTINUE 


/*  COMPUTE  THE  DEGREE  OF  FREEDOM  FOR  CHI -SQUARE 

DF  =  0 

DO  60  I  =  1,  N 

DF  =  DF  +  SAMPLE(I) 
60   CONTINUE 

DF  =  2  *  DF 


/*  COMPUTE  CHI -SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF  */ 

DO  200  LL  =  1,  2 

IF  (  DF  . EQ.  1  )  THEN 
P  =  ALPHA(LL)  /  2 
CALL  ZTABLE(P,  ZALPHA(LL)) 
CHISQR  =  ZALPHA(LL)  **  2 

ELSE  IF  (  DF  .EQ.  2  )  THEN 

CHISQR  =  -2  *  ALOG(ALPHA(LL)) 

ELSE  IF  (  DF  .  GE.  3  )  THEN 

CHISQR  =  DF  *  (1  -  2./(9  *  DF)  + 

ZALPHA(LL)  *  SQRT(  2./(9  *  DF)))  **  3 

END  IF 

LAMMAX  =  0 
LAMB I G  =  0 

DO  180  ITER  =  1,  NN 

DO  140  I  =  1,  N 

/*  GENERATE  EXPONENTIAL  RANDOM  NUMBER  WITH  MU=1  */ 
CALL  LEXPN(ISEED,  EXPO,  SAMPLE(I),  1,  0) 

/*  CONVERT  EXPONENTIAL  RANDOM  NUMBER  WITH  MU=1  */ 
/*  TO  MU=LAMBDA  AND  GET  WE I BULL  RANDOM  NUMBERS  */ 
DO  80  J  =  1,  SAMPLE(I) 

XX(I,J)  =  EXPO(J)  /  LAMSTR(I) 

W(I,J)   =  XX(I,J)  **  (1/BETA(I)) 

Y(I,J)   =  ALOG(W(I,J)) 

WEIB(J)  =  W(I,J) 
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B(J)    =  Y(I, J) 
80  CONTINUE 


/*  SUBROUTINE  SHSORT  WILL  SORT  ARRAY  B  IN  */ 
/*  ASCENDING  ORDER  BY  SHELL  SORT  ALGORITHM  */ 
CALL  SHSORT(B,  KEY,  SAMPLE(I)) 

/*  FINDING  BETA  HAT  (BETHAT)  BY  THE  FORMULA      */ 

SUM1  =  0 
SUM2  =  0 

DO  100  J  =  1,  S(I) 

SUM1  =  SUM1  +  B(J) 
100  CONTINUE 

DO  120  J  =  S(I)+1,  SAMPLE(I) 
SUM2  =  SUM2  +  B(J) 
120  CONTINUE 

BETHAT(I)  =  (SAMPLE(I))  *  KSUBM(I)  / 
*  (S(I)*SUM2  /  (SAMPLE(I)  -  S(I))  -  SUM1) 

TEMP  =  BETA(I) 

/*  SUBROUTINE  NEWTON  WILL  COMPUTE  THE  ESTIMATE  */ 
/*  OF  BETA  (BETNEW)  BY  NEWTON -RAPHSON  METHOD  */ 
CALL  NEWTON(WEIB,  SAMPLE(I),  TEMP,  BETNEW(I)) 

XTOBET(I)  =  X(I)  **  BETHAT(I) 
XBENEW(I)  =  X(I)  **  BETNEW(I) 

XHAT(I)  =  0. 
XNEW(I)  =  0. 

DO  130  J  =  1,  SAMPLE(I) 

XIJHAT(I,J)  =  W(I,J)  **  BETHAT(I) 
XHAT(I)     =  XHAT(I)  +  XIJHAT(I,J) 

XIJNEW(I,J)  =  W(I,J)  **  BETNEW(I) 
XNEW(I)     =  XNEW(I)  +  XIJNEW(I,J) 
130  CONTINUE 

LAMHAT(I)  =  SAMPLE(I)  /  XHAT(I) 
LAMNEW(I)  =  SAMPLE(I)  /  XNEW(I) 

IF  (LAMHAT(I)  . GT.  LAMMAX  )  THEN 

LAMMAX  =  LAMHAT(I) 
END  IF 

IF  (LAMNEW(I)  .GT.  LAMBIG  )  THEN 

LAMB I G  =  LAMNEW(I) 
ENDIF 
140  CONTINUE 

SUM3   =  0 
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SUM3N  =  0 
SUM4  =  0 
SUM4N  =  0 

DO  160  I  =  1,  N 

RATIO(I)   =  LAMHAT(I)  /  LAMMAX 
RATNEW(I)  =  LAMNEW(I)  /  LAMBIG 
SUM3  =  SUM3  +  RATIO(I)  *  XHAT(I) 
SUM4  =  SUM4  +  RATIO(I)  *  XTOBET( I ) 

SUM3N  =  SUM3N  +  RATNEW(I)  *  XNEW(I) 
SUM4N  =  SUM4N  +  RATNEW(I)  *  XBENEW(I) 


160  CONTINUE 


LAMMU     =  CHISQR  /  (  2  *  SUM3  ) 
RSL(ITER)  =  EXP(  -  LAMMU  *  SUM4  ) 

LAMUNR    =  CHISQR  /  (  2  *  SUM3N  ) 
RSLNEW(ITER)  =  EXP(  -  LAMUNR  *  SUM4N  ) 


180      CONTINUE 

CALL  SHSORT(RSL,  KEY1,  NN) 
CALL  SHSORT(RSLNEW,  KEY2 ,  NN) 

KK  =  (  1  -  ALPHA(LL))  *  NN 

/*  RHAT  IS  THE  RELIABILITY  OF  THE  SERIES  SYSTEM  COMPUTED   */ 
/*  BY  THE  FORMULA  */ 

/*  RHTNEW  IS  THE  RELIABILITY  OF  THE  SERIES  SYSTEM  COMPUTED  */ 
/*  BY  THE  NEWTON  -  RAPHSON  METHOD  */ 

RHAT(CASE)    =  RSL(KK) 
RHTNEW(CASE)  =  RSLNEW(KK) 

/*  SUBROUTINE  FINDJ  WILL  FIND  THE  TRUE  CONFIDENCE  LIMIT    */ 
CALL  FINDJ(RSL,  NN,  RSYS(JJ),  Jl) 
CALL  FINDJ(RSLNEW,  NN,  RSYS(JJ),  J2) 

CLR(CASE)     =  Jl  /  FLOAT(NN) 
CLRNEW(CASE)  =  J 2  /  FLOAT(NN) 

CASE  =  CASE  +  1 

200    CONTINUE 

220  CONTINUE 

DO  240  I  =  1,  N 

WRITE(6,500)  BETA(I),  SAMPLE(I) 
240  CONTINUE 

WRITE(6,550)  N,  NN 
WRITE(6,600) 
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WRITE(6,650) 
WRITE(6,670) 

CASE  =  1 

DO  280  JJ  =  1,  2 

DO  260  LL  =  1,  2 

WRITE (6, 700)  N,RSYS( JJ) ,ALPHA(LL) ,RHAT(CASE) ,RHTNEW(CASE) , 
*  CLR(CASE),  CLRNEW(CASE) 

CASE  =  CASE  +  1 

260      CONTINUE 
280  CONTINUE 


500  FORMAT(T10,' INITIAL  BETA 

T45, 'NUMBER  OF  SAMPLE  SIZES 
550  FORMAT(T 10,' NUMBER  OF  COMPONENTS 

/,T10,' NUMBER  OF  REPLICATIONS 


' ,3X,F4. 1, 
\2X,I3) 
',2X,I3, 
',2X,I3) 


600  FORMATC ' l' ,5(/).7X ,'****  RELIABILITY  OF  SERIES  SYSTEM  *****) 
650  FORMATC ///,T35,fRl,  CL1  ;  USING  THE  FORMULAR' , 

/,T35  'R2,  CL2  ;  USING  THE  NEWTON  RAHPSON  METHOD') 
670  FORMATC/// ,T6, 'NUMBER  OF' .T19 ,  ' RELIABILITY' ,T33, 'ALPHA' ,T46, 'Rl' , 

T56. 'R2' ,T66,'TRUE* ,T74,'TRUE' ,/, 

T6, 'COMPONENTS'  ,T19,' OF  SYSTEM'  ,T66,'C.  L.  l',T74,'C.L.  2', 

/,T5,76('-')) 
700  FORMATC /,T8, 15 ,T22,F5.  3 ,T33 ,F5.  2 ,T43,F7.  4 ,T53 ,F7.  4,T65 , 

F7.4,T73,F7.4) 


STOP 
END 


SUBROUTINE  ZTABLE( ALPHA, Z ALPHA) 
C 

C     «  SUBROUTINE  ZTABL2  COMPUTES  RIGHT  PERCENT  POINT  ZALPHA  FROM  » 
C     «  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  » 

C 

REAL  ALPHA,  ZALPHA 


IF  (C  ALPHA  .  GT.  0.0)  .OR.  (ALPHA  .  LT.  1))  THEN 

W  =  -  AL0GC4  *  ALPHA  *  (  1  -  ALPHA  )) 

ZALPHA  =  SQRTCW  *  (2.06118  -  (  5.72622  /  (W  +  11.6406)))) 

IF  (ALPHA  .GT.  0.5)  THEN 

ZALPHA  =  -  ZALPHA 
END  IF 
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ENDIF 

RETURN 
END 


ft  * 

*  Vr  ?V  VrV r  Vr  iV  Vr  ?V  VoWr  Vw  r  V«V  *  Vr  ^V  vV  ?V  Vc  Vr  V^ 

SUBROUTINE    SHSORT( A,KEY,N) 
DIMENSION  A(N),KEY(N) 
Ml=l 
6   M1=M1*2 

IF   (Ml    .LE.    N)      GO  TO  6 
Ml=Ml/2-l 
MM=MAX0( Ml/2,1) 
GO  TO   21 

20  MM=MM/2 

IF    (MM   .LE.     0)      GO  TO    100 

21  K=N-MM 

22  DO    1   J=1,K 
II=J 

11    IM=II+MM 

IF   (A(IM)    .GE.    A(II))      GO  TO    1 
TEMP=A(II) 
IT=KEY(II) 
A(II)=A(IM) 
KEY(II)=KEY(IM) 
A(IM)=TEMP 
KEY(IM)=IT 
I 1=1 I -MM 

IF    (II    . GT.     0)      GO  TO    11 
1    CONTINUE 
GO   TO   20 
100   RETURN 
END 


ju  ju  ju  ju  ju  ju  ju  j-  ju  ju  ju  ju  J-  -*-  ju  j-  j-  j-  ju  j^  j-  ju  j-  *V  Vr  Vr  ~V  */r  Vr  *&  *V  Vc  V"  Vc  V"  * V  *V  Vc  *V  Vr  Vr  Vr  Jr  *V  V*  ~V  *V  Vc  ~V  ~V  *V  Vr  Vr  *V  *V  V  r  "V  *  c  *  V  *'r  ~V  Vr  *V  *V  V*  *V  iV  ~V  *V  Vr  V*  Vr 

»v  * 

ju  ju  ju  j*  ju  ju  ju  ju  ju  y-  Ju  ju  ju  ju  ju  ju  y-  ju  ju  ju  ju  j-  j-  Jf  ju  j-  j-  j-  ju  j-  j-  j-  ju  j-  j-  j-  j-  j*  ju  ju  j.  j-  j-  ju  j,  ju  j-  JU  ju  ju  j-  ju  ju  j.  j-  J.  j.  ju  ju  ju  y-  ju  ju  j.  ju  ju  ju  ju  ju  j.  j-  j, 

SUBROUTINE   NEWT0N(TIME,    M,    BETA,    BETNEW    ) 
C 

C  «   SUBROUTINE  NEWTON  WILL  COMPUTE  THE  ESTIMATE   OF  BETA   (BETHAT)    » 

C  «   BY  NEWTON    -   RAPHSON  METHOD  » 

C 

REAL  TIME(M),    LNTIME(IOO),    A(3),    LNTS 

INTEGER        R 


DATA  TS    /    6    / 

C  =0 

ITER  =   0 
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R         =  M 

LNTS   =   ALOG(TS) 

DO   20    I   =    1,    R 

LNTIME(I)    =  ALOG(TIME(I)) 
C   =  C   +   LNTIME(I) 
20   CONTINUE 

C   =   C    /   R 

30   DO   60   J  =    1,    3 
SUM  =   0 

DO   40   K  =    1,    R 

IF   (    J   .EQ.     1    )    THEN 

SUM  =   SUM  +  TIME(K)    **   BETA 
ELSE 

SUM  =   SUM  +   ((TIME(K)    **   BETA)*(LNTIME(K)    **    (J-l))) 

END  IF 
40  CONTINUE 

A(J)    =   SUM 

60   CONTINUE 

/*   FUNCTION  FPRIME    IS   THE   DERIVATIVES   OF   FUNCTION  F  */ 

QUOT  =  A(2)    /   A(l) 

FPRIME   =  A(3)    /   A(l)    -    (QUOT**2)    +   ( ( 1/BETA)**2) 

F  =  QUOT   -    (1   /   BETA)    -   C 

/-•    BETA   IS   UPDATED  EACH  TIME  AND  CHECK   IF   IT  CONVERGES  */ 

BETA   =   BETA    -    F   /   FPRIME 

ITER  =   ITER  +   1 

IF    (BETA   .  GT.     25.  )    GOTO    100 

IF   (    ABS(F)    .  GT.    0.0001    )    GOTO   30 

ALPHA     =   (A(l)    /   R)   **   (1/BETA) 
BETNEW  =   BETA 

RETURN 

100  WRITE(6,*)      'DID  NOT  CONVERGE* 

WRITE (6,*)       'TRY  AGAIN  WITH  BETTER  ESTIMATE   OF   BETA' 

RETURN 
END 


yf  «»-  ..»-  j  -  y-  -'-  y  -  jTf  -y  -'-  j^  y*  ju  jj.  Jm  j,  *'r  .j,  y-  -'.  y-  j,  j-  y-  y-  j-  .'-  -'-  J  -  -'-  -y  -y  J-  yf  j-  J-  -'-  J-  j-  -y  y-  y-  yf  y,  jj.  -'..  jf  .u  -l.  y-  yf  -T-  y-  ..t-  -f-  j*  -i-  yf  -'..  -•-  y  .  y-  y-  y-  .'-  y-  y-  y  -  y-  -t-  «tf  »»f 
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SUBROUTINE  FINDJ(A,  NN,  R,  J) 
C 

C     «  SUBROUTINE  FIND J  FINDS  THE  INDEX  OF  ARRAY  A  WHICH  THE      » 
C     «  VALUE  OF  IT  IS  CLOSEST  TO  R.  » 

C 

REAL  A(NN),  R,  VALUE 

INTEGER  J 

VALUE  =  ABS(A(NN)  -  R) 
DO  100  I  =  NN-1,  1,  -1 

IF  (ABS(A(I)  -  R)  . LT.  VALUE)  THEN 

VALUE  =  ABS(A(I)  -  R) 
ELSE 

J  =  I  +  1 

RETURN 
ENDIF 

100  CONTINUE 

RETURN 
END 
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APPENDIX  D.     FORTRAN  CODE  FOR  INTERVAL  ESTIMATION 
PROCEDURE  -  NORMAL  CASE 

PROGRAM  NORMAL 
C 

C       *  * 

C  *  THIS  PROGRAM  DETERMINE  THE  ACCURACY  OF  AN  APPROXIMATE 

C  *  LOWER  CONFIDENCE  BOUND  FOR   P(  X  >  Y  ).  * 

C 

C  *  PROGRAM  NORMAL  IS  THE  CASE  WHEN  Y  IS  GIVEN  A  VALUE  YO 

C  &  X  IS  NORMALLY  DISTRIBUTED  WITH 

C  *  UNKNOWN  PARAMETERS.  * 

C  *  * 

c 

c 

REAL  YO,  P,  ZP,  MUX,  SUMX,  SUMX2 ,  XBAR 

REAL  TEMPI,  TEMP2,  TEMP3,  SIGHAT 

REAL  R(3),  SIGMAX(3),  ALPHA(2),  X(100),  Y( 100) 

REAL  RL1(54,1000),  ARLl(lOOO),  BRLi(54),  KEY(IOOO),  CIRL1(54) 

REAL  T1RL1(54),  T1CI1(54),  T2RL1(54),  T2CI1(54) 

REAL  TRLK2),  TCI1(2) 

REAL  ZHAT(IOOO),  ZKNIFE(54),  ZVAR(54),  ZBAR(54) 

REAL  MUJUNK(54) 

INTEGER  NUMX(3),  NUMY(3),  CASE,  LINE 


DATA  R  /. 95,. 99,. 995/,   SIGMAX  /. 5  , 1 ,  20/ ,   NUMX  /10,25,75/ 

DATA  NUMY  / 10, 25, 75/  ,   ALPHA  /  .2,  .1/ 

DATA  ISEED,  JSEED  /  4875,7981  / 

DATA  YO  /  400  / 

CASE  =  0 

/*  II  IS  THE  INDEX  FOR  R.  R(l)  =  .95,  R(2)  =  .99,  R(3)  =  .995   */ 
DO  500  II  =  1,  3 


C       /*  P  IS  THE  RIGTH(UPPER)  CUMULATIVE  PROBABILITY  */ 

P  =  1  -  R(II) 


C       /*  SUBROUTINE  ZTABL2  COMPUTES  RIGHT  PERCENT  POINT  ZP  */ 

C       /*  FROM  RIGHT  CUMULATIVE  PROBABILITY  P  */ 


CALL  ZTABL2(  P,  ZP  ) 


C       /*  JJ  IS  THE  INDEX  FOR  SIGMA  OF  X.  */ 

C       /--'•-   SIGMAX(l)  =  .5,  SIGMAX(2)  =  1,  SIGMAX(3)  =  3  */ 
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DO  400  JJ  =  1,  3 

MUX  =  ZP  *  SIGMAX(JJ)  +  YO 

C         /*  KK  IS  THE  INDEX  OF  NUMBER  OF  X.  */ 

C         /*  NUMX(l)  =  10,  NUMX(2)  =  25,  NUMX(3)  =  75  */ 

DO  300  KK  =  1,  3 

C  /*   LL  IS  THE  INDEX  FOR  ALPHA.  */ 

C  /*  ALPHA(l)  =  .2,  ALPHA(2)  =  .1  */ 

DO  200  LL  =  1,  2 

CALL  ZTABL2(ALPHA(LL),  ZALPHA) 

CASE   =  CASE  +  1 

TEMP  =  (NUMX(KK)  -  1.  )  /  NUMX(KK) 

TEMPO  =  SQRT(TEMP) 

MUJUNK(CASE)  =  MUX 

C  /*  REPLICATE  1000  TIMES  FOR  EACH  CASES  */ 

DO  100  1=1,  1000 

C  /*  USE  NORMAL  RANDOM  NUMBER  GENERATOR  TO  GET  */ 

C  /*  NUMX(KK)  NUMBER  OF  X.  */ 

CALL  LNORM( I SEED , X , NUMX( KK ) , 2 , 0 ) 

C  /*  THIS  PART  IS  TO  GET  SAMPLE  MEAN(XBAR)  AND         */ 

C  /*  SAMPLE  VARIANCE(XVAR)  */ 

DO  50  MM  =  1,  NUMX(KK) 

X(MM)   =  SIGMAX(JJ)  *  X(MM)  +  MUX 
50  CONTINUE 


C  /*  SUBROUTINE  VAR  WILL  COMPUTE  SAMPLE  MEAN  AND        */ 

C  /*  SAMPLE  VARIANCE        */ 

CALL  VAR(X,  NUMX(KK),  XBAR,  XVAR) 

C  /*  NOW  WE  COMPUTE   THE  LOWER  CONFIDENCE  BOUND         */ 

TEMPI  =  (XBAR  -  YO)  /  SQRT(XVAR  *  TEMPO) 
TEMP2  =  1.  /  NUMX(KK)  + 

((XBAR  -  YO)  **  2)  /  (2*(NUMX(KK)+1)*  XVAR) 
TEMPI  =  TEMPI  -  ZALPHA  *  SQRT(TEMP2)  *  TEMP 

C  /*  SUBROUTINE  ZTABL1  COMPUTE  RIGHT  CUMULATIVE  PROBA-   */ 

C  /*  BILITY  FROM  RIGHT  PERCENT  POINT.  */ 

CALL  ZTABL1(TEMP1,  ARL1(I)) 
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ARLl(I)  =  1.  -  ARL1(I) 

CALL  JKNIFE(X,  NUMX(KK),  YO ,  ZHAT(I)) 

100         CONTINUE 

CALL  VAR(ZHAT,  1000,  ZBAR(CASE),  ZVAR(CASE)) 

SIGHAT  =  l./NUMX(KK)  +  ZBAR(CASE)**2  /  ( 2*(NUMX(KK)+1) ) 

SI GHAT  =  SQRT( SIGHAT)  *  TEMP 

SIGHAT  =  ZBAR(CASE)  /  TEMPO  -  SIGHAT  *  ZALPHA 

CALL  ZTABL1( SIGHAT,  ZKNIFE(CASE) ) 

ZKNIFE(CASE)  =  1  -  ZKNIFE(CASE) 

C  /*  NON-IMSL  LIBRARY  'SHSORT'  WILL  SORT  ARL1  BY  SHELL    */ 

C  /*  SORT  ALGORITHM.  */ 

CALL  SHS0RT(ARL1,  KEY,  1000) 

DO  150  I  =  1,  1000 

RL1(CASE,I)  =  ARLl(I) 
150         CONTINUE 


C  /*  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  ARL20  WHICH  THE  */ 

C  /*  VALUE  OF  IT  IS  CLOSEST  TO  R.  */ 

CALL  FINDJ(ARL1,  R(II),  J) 

C  /*  DIVIDING  THE  INDEX  BY  1000  WILL  GIVE  US  TRUE         */ 

C  /*  CONFIDENCE  LEVEL.  */ 

CIRLl(CASE)  =  J  /  1000. 

MM  =  1000  *  (1  -  ALPHA(LL)) 
BRLl(CASE)  =  RL1(CASE,MM) 

CALL  TNYVAL(R(II),  SIGMAX(JJ),  NUMX(KK) , 
*  ALPHA (LL),  ZALPHA,  MUX,  YO ,  TRL1,  TCI1) 

TlRLl(CASE)  =  TRLl(l) 

TlCIl(CASE)  =  TCI1(1) 

T2RL1(CASE)  =  TRL1(2) 

T2CIKCASE)  =  TCI1(2) 


200       CONTINUE 
300     CONTINUE 

400    CONTINUE 
500  CONTINUE 
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WRITE(6,600) 
WRITE(6,630) 
I  =  1 
LINE  =  0 

DO  550  II  =  1,  3 
DO  550  JJ  =  1,  3 
DO  550  KK  =  1,  3 

IF  (LINE  . GE.  3)  THEN 

WRITE(6,800) 

LINE  =  0 
END  IF 

WRITE(6,710)  R(II),  SIGMAX(JJ),  NUMX(KK), 
*  BRLl(I),  BRLKI  +  1),  CIRLl(I),  CIRL1(I+1), 

ZKNIFE(I),  ZKNIFE(I+1), 
MUJUNK(I),  MUJUNK(I+1) 
1  =  1  +  2 
LINE  =  LINE  +  1 
550  CONTINUE 


WRITE(6,610) 
WRITE (6, 650) 
I  =  1 
LINE  =  0 

DO  560  II  =  1,  3 
DO  560  JJ  =  1,  3 
DO  560  KK  =  1,  3 

IF  (LINE  .GE.  3)  THEN 
WRITE(6,800) 

LINE  =  0 
END  IF 

WRITE(6,700)  R(II),  SIGMAX(JJ),  NUMX(KK), 
*  TlRLl(I),  T1RL1(I+1),  TlCIl(I),  T1CI1(I+1) 

1  =  1  +  2 
LINE  =  LINE  +  1 
560  CONTINUE 


WRITE(6,620) 
WRITE(6,650) 
I  =  1 
LINE  =  0 

DO  570  II  =  1,  3 
DO  570  JJ  =  1,  3 
DO  570  KK  =  1,  3 

IF  (LINE  .GE.  3)  THEN 
WRITE(6,800) 
LINE  =  0 

ENDIF 
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WRITE(6,700)    R(II),    SIGMAX(JJ),    NUMX(KK), 

T2RL1(I),    T2RLKI+1),    T2CI1(I),    T2CI1(I+1) 

1   =    1    +   2 


LINE   =   LINE   +    1 
5  70   CONTINUE 


600  FORMAT(ll,,5(/),7X,'****     Y   IS   GIVEN  A  VALUE  OF  YO     ****') 

610  FORMATC'l' ,5(/) ,7X, '****     Y   IS   GIVEN  A  VALUE  OF  YO     ***** , 

*  //.20X,'--   X   IS  TRUNCATED   NORMAL  WITH  PROBABILITY    .95    --') 
620   FORMATC'l' ,5(/) ,7X, *****     Y   IS   GIVEN  A  VALUE  OF  YO     ****', 

*  //,20X,'--   X    IS   TRUNCATED   NORMAL  WITH  PROBABILITY    .90    --') 


> 


630   F0RMAT(3(/),5X,'         CASE  ',6X,'      RL1( 1000*( 1-ALPHA)) 

*  4X,'TRUE   CONFIDENCE   LEVEL' , 5X, ' JACKNIFE ' , 
///,5X,'R  SX        N* ,3X, 'ALPHA  =   . 2' ,12X,'. 1*  ,11X,'. 2' , 

*  12X,'.  l',/,3X,13('-'),6X,25('-'),4X,21('-,),2X,16( *-'),/) 

650   FORMAT(3(/),5X,'         CASE  ',6X,'      RL1( 1000*( 1-ALPHA) ) 

4X,' TRUE   CONFIDENCE   LEVEL' ,///,5X,'R  SX        N',3X. 

'ALPHA   =   . 2' ,  12X. '. 1' ,11X, '. 2' , 12X, '. l' ,/,3X,13( ' -    ),6X, 
25('-'),4X,21('-f),/) 

700   FORMAT(3X,F4. 3 ,2X,F3. 1 , 2X, 12 ,4( 7X,F7. 4) ) 

710   FORMAT(3X,F4. 3 , IX, F4.  1 , 2X, 12 ,4( 7X,F7. 4) , IX, 2( 2X.F6.  4) , 2( 2X,F9.  5) ) 

800   FORMAT(/) 

RETURN 
END 
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SUBROUTINE    ZTABLK  ZALPHA  ,    ALPHA) 
C 

C  «   SUBROUTINE    ZTABL1   COMPUTE   RIGHT  CUMULATIVE   PROBABILITY  FROM      » 

C  «   RIGHT   PERCENT   POINT.  » 


C 


REAL   PI,    ALPHA,    ZALPHA 
PARAMETER   (    PI   =   3.  141592    ) 

IF    (    ZALPHA    .  EQ.     0.  0    )    THEN 

ALPHA  =   0.5 
ELSE 

'  ALPHA  =  EXP(    -2  *   ZALPHA  **   2   /   PI    ) 

ALPHA  =  SQRTO    -   ALPHA  *   (1  +  2*(  PI-3)*ZALPHA**4/(  3*PI**2) ) ) 

ALPHA  =  0. 5   *   (    1    -   ALPHA   ) 

IF   (    ZALPHA   .  LT.    0.0    )   THEN 
ALPHA  =    1    -   ALPHA 
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ENDIF 
END  IF 

RETURN 
END 


SUBROUTINE   ZTABL2( ALPHA, ZALPHA) 
C 

C  «   SUBROUTINE    ZTABL2   COMPUTES   RIGHT   PERCENT  POINT   ZALPHA  FROM        » 

C  «   RIGHT  CUMULATIVE   PROBABILITY  ALPHA  » 


C 


REAL  ALPHA,    ZALPHA 

IF   ((ALPHA   .  GT.    0.0)    .OR.     (ALPHA   .  LT.     1))   THEN 

W  =   -   AL0G(4   *  ALPHA  *   (    1    -   ALPHA    )) 

ZALPHA  =  SQRT(W  *   (2.06118    -    (    5.72622   /   (W  +   11.6406)))) 

IF   (ALPHA    .  GT.     0.  5)    THEN 

ZALPHA  =    -    ZALPHA 
ENDIF 
ENDIF 

RETURN 
END 
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SUBROUTINE   TTABLE( ALPHA,    NU ,    TALPHA) 
C 

C  «   SUBROUTINE   TTABLE   COMPUTES   RIGHT   CUMULATIVE   PROBABILITY  » 

C  «  TALPHA  FROM  T-DISTRIBUTION   BY   GIVEN   ALPHA  AND  NU  » 

C 

REAL  ALPHA,    ZALPHA,    TALPHA,    Al,    A2 ,    A3 

INTEGER  NU 

IF    ((ALPHA   .GT.     0.0)    .OR.     (ALPHA   .  LT.     1))    THEN 

CALL   ZTABL2( ALPHA,    ZALPHA) 

Al  =   (ZALPHA**2  +   1)    /   4 

A2  =   (5*ZALPHA**4  +   16*ZALPHA**2  +   3)    /   96 

A3  =   (3*ZALPHA**6  +   19*ZALPHA**4  +   17*ZALPHA**2    -    15)    /   384 

TALPHA  =  ZALPHA  *   (1  +  Al/NU  +  A2/NU**2  +  A3/NU**3) 

IF   (    ALPHA   .  GT.     0.  5    )    THEN 

TALPHA  =   -   TALPHA 
ENDIF 


4S 


END  IF 

RETURN 
END 


iV  }WoW?  Vr  VoV  Vr  V- Vr  sV  ?V  Vr  Vr  Vr  VoV  Vr  Vc  Vr  iV  ?V  V«^^ 

*  * 

iV  sV -,Wr  *  Vr  vV  Vr -/V  Vc -,V  Vr  VrV r  VoV  Vc  VoV  Vr  ?V  Vc  V^ 

SUBROUTINE  SHSORT( A,KEY,N) 
DIMENSION  A(N),KEY(N) 
Ml=l 
6  M1=M1*2 

IF  (Ml  . LE.  N)   GO  TO  6 
Ml=Ml/2-l 
MM=MAX0( Ml/2,1) 
GO  TO  21 

20  MM=MM/2 

IF  (MM  . LE.  0)   GO  TO  100 

21  K=N-MM 

22  DO  1  J=1,K 
II=J 

11  IM=II+MM 

IF  (A(IM)  .GE.  A(II))   GO  TO  1 
TEMP=A(II) 
IT=KEY(II) 
A(II)=A(IM) 
KEY(II)=KEY(IM) 
A(IM)=TEMP 
KEY(IM)=IT 
I 1=1 I -MM 

IF  (II  . GT.  0)   GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 


SUBROUTINE  VAR(Z,  NUMZ,  ZBAR ,  ZVAR) 
C 

C     «  SUBROUTINE  VAR  WILL  COMPUTE  THE  MEAN  AND  THE  VARIANCE        » 
C     «  OF  ARRAY  Z  » 

C 

REAL  Z(NUMZ),  SUMZ ,  SUMZ2,  ZBAR,  ZVAR 

SUMZ   =  0 
SUMZ 2  =  0 

DO  100  I  =  1,  NUMZ 

SUMZ   =  SUMZ   +  Z(I) 
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SUMZ2   =   SUMZ2   +   Z(I)    **   2 

100   CONTINUE 

ZBAR  =   SUMZ      /   NUMZ 

ZVAR  =  ABS(    SUMZ2    /   NUMZ    -    ZBAR  **   2    ) 

RETURN 
END 
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SUBROUTINE  FINDJ(A,    R,    J) 
C 

C  «   SUBROUTINE   FINDJ  FINDS   THE    INDEX   OF  ARRAY  A     WHICH  THE  » 

C  «  VALUE   OF   IT   IS   CLOSEST  TO  R.  » 

C 

REAL  AC  1000),    R,    VALUE 

INTEGER   J 

VALUE   =  ABS(A(1000)    -   R) 
DO    100    I   =   999,    1,    -1 

IF    (ABS(A(I)    -   R)    . LT.    VALUE)    THEN 

VALUE  =  ABS(A(I)    -   R) 
ELSE 

J  =   I   +   1 

RETURN 
ENDIF 

100   CONTINUE 

RE-1"  RN 
END 
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SUBROUTINE   TNYVAL( R , SIGMAX , NUMX , ALPHA , ZALPHA , MUX , YO , 

TBRL1,TCIR1) 
C 

C  «   SUBROUTINE   TNYVAL   IS   THE   CASE  WHEN  X   IS   TRUNCATED  NORMAL  » 

C  «   AND   Y   IS   GIVEN  A  VALUE  » 

C 

REAL  R,    SIGMAX,    ALPHA,   TALPHA,    MUX,    P(2),    ZP(2),    A(2) 

REAL  X(100),    XI,    TEMPI,    TEMP2 ,    YO 

REAL  TARLl(lOOO),    KEY(IOOO),    TRLl(lOOO),    TBRL1(2),    TCIR1(2) 

INTEGER  NUMX,    CASE,    COUNT 
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DATA  P  /  .  95,  .  90  / 

DATA  ISEED,  JSEED  /  4875,  7981  / 

DO  400  I  =  1,  2 

P(I)  =  1  -  P(I) 
COUNT  =  0 

CALL  ZTABL2(P(I),  ZP(I)) 

A(I)  =  MUX  +  ZP(I)  *  SIGMAX 

DO  200  J  =  1,  1000 

100       CALL  LNORM( ISEED,  XI,  1,  2,  0) 

XTEMP  =  MUX  +  SIGMAX  *  XI 

IF  (XTEMP  .LE.  A(I))  THEN 

COUNT    =  COUNT  +  1 

X( COUNT)  =  XTEMP 
ELSE 

GOTO  100 
END  IF 

IF  (COUNT  . LT.  NUMX)  GOTO  100 

CALL  VAR(X,  NUMX,  XBAR,  XVAR) 

C  /*  NOW  WE  COMPUTE  THE  LOWER  CONFIDENCE  BOUND  */ 

TEMPI  =  (XBAR  -  YO)  /  SQRT(XVAR*(NUMX-1)  /  NUMX) 
TEMP2  =  1.  /  NUMX  + 

((XBAR  -  YO)  **  2)  /  (2*(NUMX+1)  *  XVAR) 
TEMPI  =  TEMPI  -  ZALPHA  *  SQRT(TEMP2) 

C  /*  SUBROUTINE  ZTABL1  COMPUTE  RIGHT  CUMULATIVE  PROBA-      */ 

C  /*  BILITY  FROM  RIGHT  PERCENT  POINT.  */ 

CALL  ZTABLKTEMP1,  TARLl(J)) 

TARLl(J)  =  1.  -  TARLl(J) 

200    CONTINUE 

C        /*  NON-IMSL  LIBRARY  'SHSORT*  WILL  SORT  TARL1  BY  SHELL        */ 
C        /*  SORT  ALGORITHM.  */ 

CALL  SHS0RTCTARL1,  KEY,  1000) 

DO  300  J  =  1,  1000 
TRLl(J)  =  TARLl(J) 
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300  CONTINUE 


C  /*   SUBROUTINE   FIND J  FINDS   THE    INDEX   OF  TARL1   WHICH  THE  */ 

C  /*   VALUE   OF    IT   IS   CLOSEST  TO   R.  */ 

CALL  FINDJ(TARL1,    R,    J) 

C  /*   DIVIDING  THE    INDEX   BY    1000  WILL  GIVE   US   TRUE  */ 

C  /••-   CONFIDENCE   LEVEL.  */ 

TCIRl(I)    =  J   /    1000. 

MM  =   1000  *   (1    -   ALPHA) 
TBRLl(I)    =  TRLl(MM) 

400   CONTINUE 

RETURN 
END 
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SUBROUTINE   JKNIFE(X,    N,    YO,    ZHAT) 
C 

C  «   THIS   ROUTINE    IS   FOR  THE    'JACKNIFE'    METHOD.  » 

C  «   DUMMY   PARAMETER    IS    ZHAT.  » 


C 


REAL  X(N) ,XOMIT( 100) ,XOMBAR( 100) ,XOMVAR( 100) ,ZOMHAT( 100) , ZHAT 

INTEGER   M,    N 

M  =  N    -    1 
SUM  =   0. 

DO   200    I    =    1,    N 

DO    100   J  =    1,    M 

IF   (J   .  GE.     I)    THEN 

XOMIT(J)    =  X(J+1) 
ELSE 

XOMIT(J)    =  X(J) 
END  IF 

100  CONTINUE 

CALL  VAR(XOMIT,    M,    XOMBAR( I ) ,    XOMVAR(I)) 

ZOMHAT(I)    =   (    XOMBAR(I)    -   YO    )    /    SQRT(XOMVAR( I) ) 
SUM  =   SUM  +ZOMHAT(I) 

200   CONTINUE 

ZHAT  =   SUM   /    N 
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RETURN 
END 
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